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1 .  INTRODUCTION 


Most  diffusion  models  applied  by  che  Department  of  Defense,  the 
Environmental  Protection  Agency  (EPA),  and  other  government  agencies  are 
based  on  the  30-year-old  Pasquill  stability  category  methodology,  whose  six 
stability  classes  are  determined  as  a  function  of  wind  speed,  solar  elevation 
angle,  and  cloudiness.  Slightly  different  methods  are  used  for  rural  and 
urban  areas,  but  there  is  no  continuous  variation  between  the  two  surface  or 
between  stability  classes  for  the  same  surface.  As  a  result,  predicted 
pollutant  concentrations  are  discontinuous  functions  of  ‘he  six  stability 
classes  and  two  types  of  surfaces. 

In  1977,  the  American  Meteorological  Society  convened  a  workshop  on 
stability  classification  schemes  which  recommended  that  advances  in  boundary 
layer  theory  be  used  to  revise  the  discrete  stability  class  methodology 
(Hanna  et  al. ,  1977).  Recent  advances  in  convective  scaling  have  since 
developed  to  provide  a  better  basis  for  analyzing  dispersion  in  the 
convective  boundary  layer.  This  Phase  I  research  is  intended  to  outline  the 
structure  of  a  new  methodology  for  predicting  continuous  stability  indicators 
based  on  the  current  state-of-the-science.  The  algorithms,  intended  for 
implementation  on  a  personal  computer,  require  only  routinely  available 
meteorological  observations  and  land  use  categorizations.  The  formulation  of 
the  model  is  general.  It  is  equally  valid  over  rural  areas  or  urban 
environments  as  long  as  the  proper  model  inputs  are  specified. 

The  model  consists  of  three  components.  The  first  component  is  a  soil 
moisture  module  which  uses  precipitation  observations  and  land  use 
characteristics  to  predict  the  important  surface  moisture  parameter.  A  new 
capacitance  method  for  determining  surface  moisture  has  been  developed  for 
this  project.  Its  characteristics  have  been  compared  against  two  more  complex 
numerical  soil  moisture  models  and  it  has  been,  used  to  simulate  moisture  in 
both  rural  and  urban  environments  for  a  three  month  period  in  order  to 
illustrate  the  technique  in  an  operational  mode.  The  topic  of  moisture 
availability  is  discussed  in  Section  2.7 
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The  second  component  of  the  model  is  a  meteorological  preprocessor  which 
uses  observations  of  wind  speed  and  direction  at  a  height  of  10  m,  cloud 
cover,  surface  characteristics,  and  surface  moisture  predicted  by  the 
capacitance  model  to  determine  the  surface  sensible  and  latent  heat  fluxes, 
surface  friction  velocity,  Monin-Obukhov  length,  and  convective  velocity 
scale.  These  variables  are  important  boundary  layer  parameters  which 
determine  atmospheric  turbulence  and  dispersion  rates.  If  mixing  heights  are 
to  be  predicted,  the  model  also  requires  temperature  sounding  data.  The 
meteorological  preprocessor  is  based  on  an  energy  balance  method,  described 
in  Sections  2.1  and  2.2.  The  parameterizations  of  the  boundary  layer 
parameters  are  described  in  Section  2.4. 


The  third  component  of  the  model  is  a  series  of  similarity  relations 
which  determine  plume  dispersion  parameters,  <r  and  cr  ,  as  a  function  of  the 

y  2 

boundary  layer  parameters.  Unlike  tradition  dispersion  curves,  the 
dispersion  parameters  vary  in  a  continuous  manner  as  a  function  of  plume 
height,  atmospheric  stability,  and  surface  characteristics.  The  similarity 
relationships  are  equally  applicable  over  urban,  suburban  and  rural  areas. 
The  dispersion  component  of  the  model  is  described  in  Section  3. 


Finally,  a  summary  of  the  conclusions  and  recommendations  for  future 
Phase  II  research  are  described  in  Section  4. 
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2.  SURFACE  ENERGY  BUDGET  AND  BOUNDARY  LAYER  PARAMETERS 


A  number  of  significant  advances  have  been  made  in  recent  years  in  the 
understanding  and  characterization  of  the  structure  of  the  planetary  boundary 
layer  (PBL)  (e.g. ,  Hanna  et  al. ,  1977;  Weil,  1985;  Briggs,  1985).  The  use  of 
appropriate  boundary  layer  scaling  parameters,  which  provides  a  continuous 
representation  of  stability  and  dispersion  in  the  atmosphere,  can  improve  the 
quality  of  dispersion  predictions  (e.g.,  van  Ulden  and  Holtsiag,  1985).  The 
principal  parameters  needed  to  describe  the  boundary  layer  structure  are  the 
surface  heat  flux,  surface  momentum  flux,  and  the  boundary  layer  height. 
Several  additional  parameters,  including  the  friction  velocity,  convective 
velocity  scale,  and  the  Monin-Obukhov  length  are  derived  from  the  primary 
parameters. 

As  part  of  the  Electric  Power  Research  Institute  (EPRI)  Advanced  Plume 
project,  Hanna  et  al.  (1986)  evaluated  several  models  for  the  prediction  of 
these  boundary  layer  parameters  from  routine ly-available  meteorological 
observations.  Two  basic  methods  are  commonly  used  to  estimate  boundary  layer 
parameters  such  as  surface  heat  and  momentum  fluxes.  The  first  method  is 
referred  to  as  the  profile  method.  It  requires  at  a  minimum  observations  of 
the  wind  speed  at  one  height  and  the  temperature  difference  between  two 
heights  in  the  surface  layer,  as  well  as  knowledge  of  the  air  temperature  and 
roughness  characteristics  of  the  surface.  Monin-Obukhov  similarity  theory  is 
then  used  to  solve  the  surface  energy  fluxes  by  iteration.  The  second 
approach,  called  the  energy  budget  method,  computes  the  surface  heat  flux  by 
parameterizing  the  various  terms  of  the  surface  energy  budget  equation. 

Hanna  et  al.  (1986)  tested  and  intercompared  four  energy  budget  models 
(Holtsiag  and  van  Ulden,  1983;  Weil  and  Brower,  1983;  Berkowicz  and  Prahm, 
1982;  and  Briggs,  1982)  and  two  profile  schemes  (two-level  tower  method  and 
four-level  tower  method).  A  major  conclusion  of  the  study  was  that  the 
energy  budget  methods  were  superior  over  land  surfaces  because  of  the  high 
sensitivity  of  the  profile  method  to  small  errors  in  the  measured  temperature 
difference.  In  addition,  the  four  different  energy  budget  methods  performed 
with  similar  skill.  For  this  reason,  the  energy  budget  method  was  selected 
as  the  basis  for  estimating  the  required  urban  energy  balance  and  boundary 
layer  parameters. 
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2.1  Urban  Surface  Energy  Budget 


The  energy  balance  at  the  surface  can  be  written  as: 

Q-*Qf  *  «h  * Q.  *  %  *  “» 

where  Q,  Is  the  net  all-wave  radiation, 
is  the  anthropogenic  heat  flux, 
is  the  sensible  heat  flux, 

Q  is  the  latent  heat  flux, 
e 

is  the  ground/storage  heat  flux,  and, 

Q  is  the  net  heat  advection. 
a 

Urbanization  changes  the  properties  of  the  surface  and  atmosphere  in 
several  important  ways  which  must  be  included  in  the  parameterization  of  the 
various  terms  of  the  energy  balance  equation.  For  example,  the  aerodynamic, 
thermal,  hydrological,  and  radiative  characteristics  of  the  surface  are 
changed  as  a  result  of  urbanization  (Oke,  1978).  These  changes  have 
important  effects  on  atmospheric  stability  and  turbulence  in  the  urban 
environment.  For  example,  buildings  and  other  obstacles  to  the  wind  flow 
create  a  very  rough  surface  which  increases  the  mechanically  generated 
atmospheric  turbulence  in  the  lowest  layer  of  the  atmosphere.  Briggs  (1983) 
estimates  increases  in  the  surface  friction  velocity  of  30-50%  in  urban  areas 
relative  to  surrounding  suburban  and  rural  areas.  The  thermal  properties  of 
the  urban  surface  are  altered  by  the  presence  of  dense  building  materials 
which  efficiently  store  heat.  In  addition,  anthropogenic  heat  sources  cam 
contribute  significantly  to  the  urban  energy  budget.  The  large  fraction  of 
the  urban  surface  impervious  to  water  and  urban  runoff  routing  systems  affect 
the  available  moisture  and  therefore,  the  partitioning  of  available  energy 
into  sensible  and  latent  heat  fluxes.  Godowitch  et  al.  (1981),  White  et  al. 
(1978),  Hildebrand  and  Ackerman  (1984)  and  others  report  urban  sensible  heat 
fluxes  2-4  times  higher  than  those  in  nearby  rural  areas.  In  addition,  urban 
areas  tend  to  have  higher 'concentrations  of  atmospheric  pollutants,  which 
alters  the  spectrum  and  intensity  of  solar  radiation  received  at  the  surface 
(Landsberg,  1981). 
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Cleugh  and  Oke  (1986)  present  simultaneous  energy  balance  observations  at 
rural  and  suburban  sites  in  Vancouver.  B. C.  which  illustrate  some  of  the 
effects  of  urbanization.  The  major  components  of  the  energy  balance  equation 
are  shown  in  Figure  1  averaged  over  a  period  of  30  days  during  the  summer  of 
1983  for  both  sites.  In  the  suburban  area,  the  net  radiation  increased  4%, 
the  sensible  heat  flux  increased  51%,  and  the  latent  heat  flux  decreased  46% 
as  compared  to  the  rural  values.  The  Bowen  ratio  (ratio  of  sensible  to 
latent  heat  flux)  averaged  0.5  in  the  rural  area  and  1.3  in  the  suburban 
area,  reflecting  decreased  moisture  availability  in  the  suburban  environment. 


The  energy  budget  methodology  provides  a  general  framework  for 
predicting  energy  fluxes,  and  profiles  of  winds,  temperature,  and  turbulence 
which  is  valid  for  all  types  of  underlying  surfaces.  By  properly  defining 
the  surface  characteristics  and  related  parameters,  the  energy  budget 
technique  can  be  used  to  derive  estimates  of  atmospheric  stability  and  plume 
dispersion  parameters  which  are  continuous  in  both  space  and  time  for 
surfaces  ranging  from  highly  urbanized  to  suburban  to  rural  land  use. 

2.2  Parameterization  of  the  Terms  of  the  Surface  Energy  Budget  Equation 

The  solution  of  the  energy  budget  equation  is  based  primarily  on  the 
methods  of  Holtslag  and  vein  Ulden  (Holtslag  and  van  Ulden,  1983;  van  Ulden 
and  Holtslag,  1985).  The  energy  budget  equation  is  solved  for  the  sensible 
heat  flux  by  parameterizing  the  other  terms  of  the  equation  in  terms  of  known 
quantities.  Meteorological  preprocessors  implementing  these  techniques  for 
urban  areas  have  been  developed  and  applied  by  Scire  et  al.  (1986)  (used 
predicting  plume  dispersion  and  dry  deposition  rates  in  New  York  City)  and 
Hanna  and  Chang  (1991  and  1992)  (as  part  of  the  HPDM-Urban  model  development 
and  evaluation  program). 

The  parameterization  of  the  various  energy  budget  terms  is  described 
below. 
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ENERGY  Fi-UX  PENSIJY  |Wm'J) 


Figure  1 


Variation  of  average  energy  balance  components  for  (a)  rural  and 

(b)  suburban  sites  in  Vancouver,  B.  C.  for  30  summer  days  in  1983 

Q*  is  net  radiation  flux,  Q,  is  sensible  heat  flux,  Q  is  latent 

n  e 

heat  flux,  and  is  ground  heat  flux.  (From  Cleugh  and  Qke, 
1986.  ) 


NET  RADIATION 


The  net  radiation,  Q#,  is  computed  as  the  residual  of  incominp 
(short-wave  plus  long-wave)  radiation  and  outgoing  (long-wave)  radiation, 
can  be  expressed  (Hoitslag  and  van  Ulden,  1983;  Landsberg,  1981)  as: 


(l-A)Q  ♦  Q.  ,  -  Q. 

sw  lw-d  lw-u 


(2) 


2 

where  is  the  incoming  short-wave  radiation  (W/m  ),  consisting  c*  a  direct 

solar  radiation  term  (Q  )  plus  a  diffuse  radiation  term  (Q  , ), 

sw-s  r  sw-d 

A  is  the  albedo  of  the  surface, 

2 

C*lw-d  is  the  incoming  long-wave  atmospheric  radiation  (W/m  ),  and, 

CL  is  the  long-wave  radiation  (W/m^)  emitted  by  the  surface, 
lw— u 


Both  components  of  the  short-wave  radiation  term  are  altered  by 

conditions  in  the  urban  environment.  Direct  solar  radiation,  Q  ,  is 

sw-s 

substantially  decreased  due  to  attenuation  by  urban  air  pollutants. 

Typically,  Qsw_s  is  about  40%  less  in  urban  areas  than  in  surrounding  rural 

areas  (Unsworth  and  Monteith,  1972).  The  ultraviolet  (UV)  portion  of  the 

spectrum  is  affected  much  more  than  the  longer  wavelengths.  However,  diffuse 

radiation,  Q  ,,  is  increased  in  urban  areas.  The  net  effect  is  that  Q  is 
sw-d  sw 

typically  10-20%  less  in  urban  areas  (Oke,  1978;  Landsberg,  1981).  Another 
mitigating  factor  in  the  net  radiation  budget  is  that  urban  areas  generally 
have  lower  surface  albedos  than  most  rural  environments,  and  thus  lower 
losses  of  short-wave  radiation  due  to  reflection  from  the  surface.  For 
example,  typical  urban  albedos  in  St.  Louis  are  0. 12-0. 13  whereas  albedos  In 
the  surrounding  rural  land  surfaces  are  0.15-0.17  (White  et  al. ,  1978). 


The  individual  long-wave  terms  of  the  net  radiation  equation  are  also 
affected  in  an  urban  area.  The  warmer  urban  surface  temperatures  increase 
the  emitted  long-wave  radiation  (Q^w  L  However,  the  incoming  long-wave 
atmospheric  radiation  is  also  increased  due  to  higher  pollutant 

levels.  The  net  effect  during  the  day  is  usually  a  slightly  smaller 
long-wave  loss  in  the  urban  area.  At  night,  the  long-wa''e  loss  is  slightly 
larger  in  the  urban  area  compared  to  rural  areas  (Oke,  1978). 


Although  each  component  of  the  net  radiation  (Q.)  equation  is  altered  by 
urban  conditions,  the  overall  net  effect  is  that  the  urban  Q.  is  usually  not 
substantially  different  (i.e. ,  is  within  about  5%)  in  mid-latitude  cities 
from  that  in  the  surrounding  rural  areas  (Oke,  1978;  White  et  al. ,  1978). 
However,  Oke  (1978)  points  out  exceptions  such  as  in  snow  covered  conditions 
or  in  cities  in  desert  areas,  where  the  rural  albedo  may  be  signif icantly 
higher.  For  example,  typical  albedos  are  in  the  range  0.  40  to  0.  95  over 
rural  snow  covered  surfaces  (Oke,  1978)  aud  0.20  to  0.45  in  desert  areas 
(Hanna  and  Chang,  1991).  Snow  in  urban  areas  is  usually  quickly  removed  by 
plowing  or  salting  operations,  and  tends  not  to  stick  to  the  vertical  faces 
of  buildings.  Also,  the  snow  melts  faster  due  to  the  urban  heat  island 
effect.  The  urban  snow  that  does  remain  tends  to  have  a  lower  albedo  due  to 
soiling  by  urban  pollution.  In  situations  such  as  these,  where  the  urban 
albedo  is  much  lower  than  that  in  the  surrounding  rural  areas,  the  term 
(l-AiQ^  can  be  substantially  higher  in  the  urban  area,  leading  to  higher 
urban  values  of  net  radiation. 


Holtslag  and  van  Ulden  (1983)  provide  the  following  parameterizations  of 
the  short-wave  and  long-wave  radiation  terms  in  the  net  radiation  equation: 

(l-A)Q  ♦  c.T6  +  c,N  -  <rT4 
sw  1  2 

Q«  =  -  (3) 

1  +  C3 

Qsw  =  (aisin*  +  a2Xl  +  b^2)  (45 


c 


3 


0.38 


r  (i-cc) cs)  ♦  i 


]  S  +  1 


(5) 


where  T 
<r 
N 

a 

S 

A 


is  the  measured  air  temperature  (deg.  K), 

-8  2  4 

Is  the  Stefrn-Boltzmann  constant  (5.67  x  10  W/m  /deg.  K  ), 
is  the  fraction  of  the  sky  covered  by  clouds, 
is  the  solar  elevation  angle  (deg. ), 
is  an  empirical  surface  moisture  parameter,  and, 

's  the  slope  of  the  saturation  enthalpy  curve  [3=s/rj,  where 

s=d(qs )/3(T)  and  y=c^/L, 

is  the  latent  heat  of  water  vaporization. 
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qg  is  the  saturation  specific  humidity,  and, 

0^  is  the  specific  heat  at  constant  pressure. 

The  four  terms  in  the  numerator  of  Eqn.  (3)  account  for  absorption  of 
short-wave  radiation  at  the  surface,  incoming  long-wave  radiation  from 
gaseous  components  of  the  atmosphere  (e.g.,  water  vapor  and  carbon 
dioxide),  incoming  long-wave  radiation  due  to  clouds,  and  outgoing  long-wave 
radiation  from  the  surface,  respectively.  The  factor  in  the  denominator 
(l+c^),  results  fr  a  the  use  of  air  temperature  rather  than  the  more 
dlff icult-to-determine  surface  radiation  temperature  in  the  equation.  The 
term  in  the  first  set  of  parentheses  in  Eqn.  (4)  represents  short-wave  solar 
radiation  in  the  absence  of  clouds.  The  second  term  ( l+b^h<  ) ,  accounts  for 
the  reduction  of  incoming  solar  radiation  due  to  clouds  ( b ^  is  negative).  The 
values  for  the  empirical  constants  c^,  c^,  a^,  a^.  b^,  and  b^  suggested  by 
Holtslag  and  van  Ulden  (1983)  are  shown  in  Table  l.  It  must  be  stressed  that 
these  "constants"  were  derived  from  observations  in  the  Netherlands  and  may  be 
expected  to  vary  somewhat  at  other  geographic  locations  characterized  by 
different  land  use  patterns,  cloud  types,  and  solar  elevation  angles. 
Furthermore,  note  that,  in  this  system,  clouds  are  parameterized  by  a  single 
number,  N,  which  is  the  fractional  cloud  cover.  Ihe  system  does  not  account 
for  type,  elevation,  and  density  of  clouds.  Perhaps  in  the  future,  field 
data  can  be  used  to  refine  these  parameter izat ions. 

The  parameter,  S,  describes  the  slope  of  the  curve  in  the  so-called 
Clausius-Clapeyron  diagram,  which  is  a  fixture  in  introductory  dynamic 
meteorology  textbooks.  Table  2  contains  values  of  S  as  a  function  of  T  (in 
5°C  increments),  for  a  pressure  of  1000  mb. 

The  term  (l+c^)  1  of  Eqn.  (3),  which  scales  the  net  radiation,  depends  on 
both  the  ambient  temperature  and  the  surface  moisture  parameter.  As  shown 
in  Figure  2,  its  sensitivity  to  soil  moisture  increases  at  higher 
temperatures.  At  35°C,  (l+c^)  *  increases  about  47%  from  very  dry  to  very  wet 
conditions  (i.e.  a  ranging  from  0.0  to  1.4),  but  increases  only  about  36% 
for  the  same  range  of  a  at  20°C,  and  19%  at  0°C.  Over  a  typical  range  of 
conditions  expected  at  a  single  site,  the  sensitivity  of  net  radiation  to 
variations  in  (1+c^)  1  is  more  of  the  order  of  20-30%. 


G 


Table  1 


Values  of  Net  Radiation  Constants 
(Holtslag  and  van  Ulden,  1983) 


Constant 


Value 


990  W/m 
-30  W/m2 

-0.  75 
3.  4 


5.31X10-13  W/m2/deg.  K6 
60  W/m2 
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Table  2 

Slope  of  the  Saturation  Enthalpy  Curve, 

S  *  £3q  /3T) (L/c  ) 
s  p 

as  a  Function  of  Temperature  for  a 
Standard  Pressure  of  1000  mb. 

(Holtslag  and  van  Ulden,  1983) 


T  (°C)  S 


-5 

0.  498 

0 

0.694 

5 

0.943 

10 

1.27 

15 

1.67 

20 

2.22 

25 

2.86 

30 

3.70 

35 

4.76 

11 


1/ (1+C3) 


Figure  2.  Plot  of  the  term  (l+c^)  *  of  the  net  radiation  equation  as  a 

function  of  the  surface  moisture  parameter  for  three  values  of 
temperature. 
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The  short-wave  radiation,  Qsw>  as  predicted  by  Eqn.  (4)  is  plotted  as  a 
function  of  time  of  day  in  Figure  3.  The  short-wave  radiation  data  is 
computed  for  a  latitude  of  40°  N  for  total  cloud  covers  of  0,  0.50.  0.75,  and 
1.0  on  two  days  of  the  year  (winter  solstice  and  summer  solstice).  Less 
than  a  4%  reduction  in  solar  radiation  is  indicated  for  a  cloudiness  factor, 
N,  up  to  0.50.  During  overcast  conditions  (N  *  1.0),  the  short-wave 
radiation  is  reduced  to  0. 25  of  the  clear-sky  value. 

Although  the  surface  albedo,  A,  is  nearly  constant  for 
elevation  angles  (i.e. ,  greater  than  30°),  it  increases  for 
(Coulson  and  Reynolds,  1971;  Iqbal,  1983).  Hanna  and  Chang 
the  use  of  the  following  relationship  based  on  Iqbal  (1983) 
the  surface  albedo  as  a  function  of  elevation  angle. 

A  =*  A'  +  ( 1-A'  )  ea*+b  (6) 

b  =  -0.5  (1-A')2  (7) 

where  A'  is  the  albedo  for  the  sun  directly  overhead, 
is  the  solar  elevation  angle  (degrees),  and, 
a  is  an  empirical  constant  (~  0.1). 


high  solar 
lower  angles 
( 1991 )  recommend 
for  estimating 


ANTHROPOGENIC  HEAT  FLUX 


In  urban  areas  with  high  population  densities  and/or  high  per  capita 

energy  usage,  the  heat  generated  by  man-made  sources  can  signif icantly  affect 

the  overall  energy  budget.  Table  3  shows  the  average  anthropogenic  heat 

flux  (Qf)  and  net  radiation  (Q,)  from  several  urban  areas.  Anthropogenic 

heat  from  sources  such  as  space  heating/cooling,  transportation,  and 

industrial  facilities  sometimes  exceeds,  on  an  annual  basis,  the  energy  input 

from  net  radiation  (e.g. ,  in  Manhattan  and  Montreal).  In  densely  populated 

northern  cities,  seasonal  values  of  Qf  sometimes  exceed  Q„  by  a  large  margin. 

*  2 

For  example,  in  Montreal,  the  average  wintertime  Qf  is  153  W/m  ,  compared  to 

2  1 
a  value  of  Q,  of  13  W/m  . 

The  patterns  of  anthropogenic  heat  generation  vary  considerably 
depending  on  climatological  factors.  The  ratio  of  wintertime  anthropogenic 
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GSW  IH/M#*2) 


TIME  (LST) 


Figure  3b.  Short-wave  solar  radiation,  Qsw  predicted  by  Eqn.  (4)  as  a 

function  of  time  of  day.  Q  is  shown  for  a  latitude  of  40°N  and 

7  sv 

total  cloud  fractions,  N,  of  0.0,  0.5,  0.75,  and  1.0  at  the  time 
of  the  winter  solstice. 


Table  3 

Average  Anthropogenic  Heat  Flux  (Qf)  and 
Net  Radiation  (Q, )  for  Several  Urban  Areas 
(From  Oke,  1978) 


Urban  area/ 

latitude/ 

period 

Population 

(x  106) 

Population 

density 

2 

(persons/km  ) 

Per  capita 

energy  usage 
(MJxl03/yr ) 

(W/m2) 

Q. 

(W/m2) 

Manhattan  (40°N) 

annual 

1.7 

28,810 

128 

117 

93 

summer 

40 

winter 

198 

Montreal  (45°N) 

annual 

1. 1 

14, 102 

221 

99 

52 

stumer 

57 

92 

winter 

153 

13 

Budapest  (47°N) 

annual 

1.3 

11,500 

118 

43 

46 

stumer 

32 

100 

winter 

51 

-8 

Sheffield  (53°N) 

annual 

0.5 

10, 420 

58 

19 

56 

West  Berlin  (52°N) 

annual 

2.3 

9,830 

67 

21 

57 

Vancouver  ( 49°N ) 

annual 

0.6 

5,360 

112 

19 

57 

stumer 

15 

107 

winter 

23 

6 

Hong  Kong  (22°N) 

annual 

3.  9 

3,730 

34 

4 

-110 

Singapore  (1°N) 

annual 

2.  1 

3,700 

25 

3 

-110 

Los  Angeles  (34°N) 

annual 

7.0 

2,000 

331 

21 

108 

Fairbanks  (64°N) 

annual 

0.03 

810 

740 

19 

18 
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heat  input  to  the  summertime  value  is  nearly  5  in  Manhattan,  but  only  1.5  in 

Vancouver,  which  has  relatively  mild  winters.  Hanna  and  Chang  (1991) 

2  2 

computed  median  values  of  36  U/ra  (summer)  and  102  W/m  (winter)  based  on 
Oke's  (1978)  data. 

Anthropogenic  heating  also  shows  considerable  spatial  variability.  For 

example,  in  Greater  London,  McGoldrick  (1980)  found  that  Q_  varied  from 

2  f  2 
0  to  5  W/m  in  outer  areas  of  the  city,  but  exceeded  100  W/m  over  several 

square  kilometers  in  the  center  of  the  city  (Landsberg,  1981).  The  maximum 

2 

value  averaged  over  one  square  kilometer  was  234  W/m  (which  is  more  than 

2 

twice  the  average  global  solar  radiation  of  106  W/m  ).  The  spatial 

2 

variability  of  the  average  annual  anthropogenic  heat  flux  for  2  km  grid 

cells  in  Tokyo  is  shown  in  Figure  4  (Kimura  and  Takahashi,  1991).  The 

2 

anthropogenic  heat  flux  ranges  from  less  than  10  W/m  over  wide  areas  of  the 

2 

metropolitan  area  to  over  100  W/m  in  the  inner  core  of  the  city. 

Hanna  and  Chang  (1991)  point  out  the  importance  of  the  elevation  of  the 
anthropogenic  energy  release.  For  example,  the  heat  input  from  a  power  plant 
occurs  at  an  effective  height  of  several  hundred  meters,  wi.xle  home  heating 
and  vehicle  exhaust  emissions  are  released  near  the  surface.  They  note  that 
only  the  near-surface  components  of  anthropogenic  heat  flux  should  be 
considered  in  the  computation  of  the  surface  energy  balance. 

These  discussions  point  out  that  the  anthropogenic  heat  input  is 
significant,  but  that  its  magnitude  and  its  spatial  and  temporal  variability 
are  rather  uncertain.  Consequently,  it  is  difficult  to  parameterize.  For 
example,  suppose  one  wanted  to  model  Minneapolis,  which  is  not  in  Table  3. 

All  one  can  say  is  that  the  average  anthropogenic  heat  flux,  Qf,  over  the 
urban  area  is  probably  about  50  w/m  ,  with  an  uncertainty  of  plus  or  minus  a 
factor  of  two.  The  uncertainty  would  be  multiplied  by  another  factor  of  two 
or  three  if  is  estimated  at  a  specific  location,  time  of  year,  and  hour  of 
day.  We  recommend  that  local  data  on  energy-use  patterns  be  analyzed  before 
making  local  estimates  of  Q^.  for  use  in  the  surface  energy  balance  equation. 
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Symbols 


x  10“  kcal/year/4km2  W/m2 
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33 

20 

X 

6.6 

30 
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a 

13. 

50 

o 

17. 

60 

A 
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30 
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0 

30. 

100 

* 

33. 

120 

a 
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140 

A 

46. 
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■ 
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a 
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A 
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a 
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Figure  4.  Average  annual  anthropogenic  heat  flux  for  2  km  square  grid  ceils 
in  the  Tokyo  Metropolitan  area  (From  Kimura  and  Takahashi,  1991). 
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GROUND/STORAGE  HEAT  FLUX 


The  flux  of  heat  into  the  ground  or  storage  in  surface  materials,  Q  ,  is 

g 

usually  parameterized  during  the  daytime  as  a  fraction  of  the  net  radiation 
(e.g, ,  DeBruin  and  Holtslag,  1982;  Oke,  1978). 


g 


°8  °* 


(s: 


where  c  is  an  empirical  coefficient  which  depends  on  the  properties  of  the 

o 

surface.  Holtslag  and  van  Ulden  (1983)  obtained  a  value  of  c  of  0. 1  for  a 

3 

grass  covered  surface  in  the  Netherlands.  Oke  (1982)  indicates  that  typical 

ranges  for  c  are  0.05  to  0.25  in  rural  areas,  0.20  to  0.25  in  suburban 
3 

areas,  and  0.25  to  0.30  in  urban  regions  and  suggests  that  typical  values  of 

c  are  0.15,  0.22,  and  0.27  for  rural,  suburban,  and  urban  areas, 
g 

respectively.  The  larger  values  of  c  for  urban  area  reflect  the  greater 

g 

ability  of  building  materials  and  pavement  to  transport  heat  (i.e. ,  greater 

thermal  conductivity,  k)  and  the  greater  heat  capacity,  C,  of  some  urban 

urban  building  materials.  Oke  (1982)  suggests  that  the  thermal  admittance, 
1/2 

H=(kC)  or  thermal  inertia  is  the  most  relevant  measure  of  the  thermal 

responsiveness  of  a  surface.  The  greater  thermal  admittance  of  urban  areas 

is  expected  to  account  for  increased  heat  storage  in  cities,  although  direct 

measurements  of  Q  on  urban  scales  are  difficult  to  make  because  of  the 
g 

variety  of  urban  materials  and  complex  geometry  of  the  surface  elements.' 


In  some  cases,  the  heat  storage  in  the  buildings  can  exceed  the  ground 
storage.  However,  when  the  buildings  have  significant  size,  the  concept  of  a 
"surface"  energy  balance  breaks  down,  and  it  is  necessary  to  calculate  the 
energy  balance  for  a  relatively  deep  layer  extending  from  the  building  tops 
down  to  the  ground  surface.  Lacking  detailed  measurements  in  this  complex 
system,  most  modelers  assume  that  the  simple  parame ter izat ions  described 
above  are  still  vaiid  approximations.  These  assumptions  are  probably 
reasonable  during  periods  of  strong  solar  insolation,  but  will  tend  to  break 
down  at  night  and  during  cloudy  conditions,  when  the  uncertainties  in  the 
energy  fluxes  are  larger  than  the  fluxes  themselves  (Hanna  and  Chang,  1992). 
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SENSIBLE  AND  LATENT  HEAT  FLUXES 


One  of  the  most  important  changes  resulting  from  urbanization  is  the 

tendency  for  a  shift  in  the  partitioning  of  available  energy  during  dry 

conditions  from  latent  into  sensible  heat.  This  shift  results  from  a  number 

of  factors,  including  decreased  surface  moisture  availability  in  the  urban 

area  due  to  the  hydrological  characteristics  its  surface,  higher  urban 

surface  temperatures,  and  less  vegetation  in  the  urban  area.  Measurements  of 

urban  sensible  heat  fluxes  two  to  four  times  those  in  the  surrounding  rural 

areas  are  common  (e.g. ,  Godowitch  et  al. ,  1981;  White  et  al.  ,  1978; 

Hildebrand  and  Ackerman,  1984).  Figure  5  plots  observations  of  the  net 

radiation,  Q,,  and  sensible  heat  flux.  Q.  ,  for  three  sites  in  the  St.  Louis 

h 

area  which  illustrates  the  increase  in  sensible  heat  with  increased 
urbanization.  Site  105  (urban  commercial  with  little  vegetation).  Site  107 
(urban  residential  with  many  trees),  and  Site  109  (.rural  agricultural).  Also 
shown  is  the  Bowen  ratio,  B,  defined  as  the  ratio  of  sensible  to  latent  heat 
fluxes,  i. e. , 

B  =  -  '  (9) 

Qe 

Oke  (1978)  gives  the  following  typical  values  and  ranges  for  the  Bowen 
ratio: 

Urban:  8  =  1.5  (range  0.5  to  over  4.0) 

Suburban:  B  =  1.0  (range  0.25  to  2.5) 

Rural:  B  =  0. 5  (range  0.1  to  1.5) 

During  wet  conditions,  i. e. ,  immediately  following  a  rain  event,  the 
differences  between  urban  and  rural  latent  heat  fluxes  are  much  less  as 
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hour  (crn 

Figure  5.  Diurnal  variati.cn  of  the  sensible  heat  flux  (H)  at  Sites  205 
(urban  commercial),  107  (urban  residential),  and  109  (rural 
agricultural)  in  the  St.  Louis  area  (top  figure).  Also  shown  is 
the  net  radiation  at  Site  105  (top  figure)  and  the  Bowen  ratio  for 
Sites  105  and  i09.  Data  averaged  over  the' month  of  August,  1976. 
(From  Ching  et  al. ,  1978). 
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latent  heat  flaxes  dominate  due  to  evaporation  from  the  wet  surface.  This  is 
a  reason  for  the  considerable  overlap  in  the  ranges  of  Bowen  ratio  for  the 
rural,  suburban  and  urban  areas.  Also,  the  presence  of  irrigation  systems 
(e.g. ,  agricultural  irrigation  in  rural  areas  or  urban  suburban  lawn 
watering)  caui  significantly  affect  the  partitioning  of  latent  and  sensible 
heat  over  the  affected  areas. 

Although  some  micrometeorological  models  ailow  specification  of  the 
Bowen  ratio,  which  determines  the  relative  magnitudes  of  the  sensible  and 
latent  heat  fluxes  directly,  Holtslag  and  van  Ulden  (1983>  recommend  a 
different  technique  based  on  a  simplified  version  of  the  Penman-Monteith 
model.  The  Penman-Monteith  expression  for  latent  heat  flux  (neglecting 
anthropogenic  heat  flux)  is  given  by  van  Ulden  and  Holtslag  (1995)  as: 


r  S 


r  p  L 


(1  ♦  r  S) 


W.  -Qg)  - 


( 1  +  r  S)  r 


[qs(T)  “  q) 


and 


r  *  r  /(r  +  r  ) 

a  a  s 


ra  -  (T  -  To)/(0.u.) 


and  from  deBruin  and  Holtslag  (1982), 
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where  r  is  the  aerodynamic  resistance  to  water  vapor  transfer  from  the 

3l 

surface  to  the  air, 

rg  is  the  surface  resistance  to  water  transfer  from  soil  and  vegetation 
to  the  surface, 

qg  is  the  saturation  specific  humidity, 
q  is  the  actual  specific  humidity  of  the  air, 

L  is  the  latent  heat  of  vaporization, 

T, Tq  are  the  air  and  surface  temperatures,  respectively. 
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0.  is  the  temperature  scale  for  turbulence  heat  transfer,  ana 
u.  is  the  surface  friction  velocity. 

The  first  terra  of  Eqn.  (10)  accounts  for  thermodynamic  evaporation,  or 

evaporation  due  to  the  presence  of  an  external  energy  source  (i.e.,  of 

strength  Q*  -  Q  ).  The  second  terra  is  called  aerodynamic  evaporation  because 
§ 

it  accounts  for  evaporation  enhancement  due  to  the  flow  of  unsaturated  air 
over  the  surface.  Because  sorae  of  the  terms  of  the  Penman-Monte ith 
equation  are  difficult  to  evaluate  operationally,  van  Ulden  and  Holtslag 
recommend  the  use  of  a  simplified  approximation  to  the  expression  (de  Bruin 
and  Holtslag,  1982;  van  Ulden  and  Holtslag,  1985)  for  latent  heat  flux 
referred  to  as  the  modified  Priestley-Taylor  model: 
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where  a  is  an  empirical  surface  moisture  parameter,  and 
B'  is  an  empirical  coefficient. 

The  basis  for  the  Eqn.  (14)  is  that  there  tends  to  be  a  strong 

correlation  of  both  the  thermodynamic  and  aerodynamic  evaporation  rates  with 

the  equilibrium  evaporation,  Q  (eq)  (defined  as  the  evaporation  which  would 

© 

occur  with  a  wet  surface  (r  *  0)  and  saturated  air  (q  *  q  ).  Thus,  Q  (eq)  = 

s  s  e 

[S/(S+1 ) ] (Q,-  Qg).  This  correlation  is  due  to  the  fact  that  (qg-q)  and 

(Q*-Q  )  exhibit  similar  diurnal  patterns  (e.g.,  see  Figure  6).  In  Eqn.  (14), 
8 

the  aerodynamic  term  is  split  into  a  portion  correlated  with  Qg(eq)  and  an 

uncorrelated  portion.  The  unccrrelated  part  is  represented  by  the  empirical 

2 

constants  aB' .  Holtslag  and  van  Ulden  (1983)  found  that  B'~  20  W/m  .  A 
simple  scheme  for  parameterizing  the  surface  moisture  parameter  is  discussed 
in  more  detail  in  the  next  section. 

In  Eqn.  (14),  the  Q^.  flux  terra  has  been  added  to  the  right  side  of 
Holtslag  and  van  Ulden’ s  equations  to  retain  the  generality  of  the  equation 
for  the  urban  case  where  anthropogenic  heat  fluxes  may  be  significant.  The 
term  a 8'  results  in  positive  latent  heat  fluxes  during  periods  when  the 
sensible  heat  flux  is  small  but  negative  (e.g.,  around  sunrise  and  sunset), 
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Figure  6.  Diurnal  variation  of  the  water  vapor  deficit  (qs~q)  and  the 
equilibrium  evaporation  (Qg(eq)).  From  de  Bruin  and 
Holtslag  (1982). 
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which,  as  illustrated  in  Figure  7,  is  an  observed  feature  of  the  behavior  of 
sensible  and  latent  heat  fluxes  (e.g. ,  Qke,  1982;  van  Ulden  and  Holtslag, 
1985). 

The  sensible  heat  flux,  is  determined  by  Holtsiag  and  van  Ulden 
(1983)  as: 


'  (1  -  a) (S)  +  1  ' 

S  +  1 


(Q.  *  Qf  *  Qg> 


(a)  O'  ) 


2.3  Moisture  Availability 


(15) 


Surface  moisture  availability  is  one  of  the  key  factors  which  determines 
the  partitioning  of  available  energy  into  sensible  and  latent  heat,  and  thus, 
it  has  an  important  effect  on  atmospheric  stability,  turbulence  levels,  and 
dispersion  rates.  For  example,  Segal  et  al.  (1990)  found  that  a  change  in 
surface  moisture  from  dry  to  wet  conditions  resulted  in  a  one  Pasquill 
stability  class  shift  toward  more  neutral  conditions  around  mid-day,  and 
larger  stability  shifts  in  the  morning  and  afternoon.  They  found 
consistently  lower  values  of  the  ratio  of  convective  velocity  scale  to  wind 
speed  (i.e. ,  w./U),  which  is  used  as  a  measure  of  stability,  with  increases 
in  soil  moisture.  Their  study  also  found  significantly  lower  mixing  heights 
(by  factors  of  1.5  to  2.0)  over  wet  surfaces. 


The  Holtslag  and  van  Ulden  model  uses  the  surface  moisture  parameter,  a, 
(sometimes  called  the  Priestley-Taylor  parameter)  as  a  measure  of  moisture 
availability.  The  relationship  between  a  and  the  Bowen  ratio,  B,  can  be 
obtained  by  combining  Eqns.  (14)  and  (15). 
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Figure  7.  Diurnal  variation  of  the  components  of  the  surface  energy  budget 

* 

measured  at  a  rural  site  in  the  Netherlands.  Net  radiation  (Q  ). 
latent  heat  flux  (LE),  sensible  heat  flux  (H),  and  soil  heat 
flux  (G).  (From  van  Ulden  and  Holtslag,  1985). 
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Holtslag  and  van  Ulden  and  others  have  applied  the  modified 
Priestley-Taylor  approach  over  a  variety  of  surfaces.  Generally,  these 
applications  involved  the  specification  of  a  as  a  constant  over  the  length  of 
the  simulation.  Typical  values  of  a,  based  on  empirical  data  of  Holtslag  and 
van  Ulden  and  summarized  by  Hanna  and  Chang  (1991)  are: 

a  *  0.2  (arid  rural  areas) 

a  a  o.  5  (urban  areas,  some  parks,  crops  and  fields  during  mid-summer 
when  rain  has  not  fallen  for  several  days) 
a  »  0.8  (crops,  fields,  or  forest  with  sufficient  moisture) 
a  =  1.0  (normal  wet  grass  in  a  moderate  climate) 

Two  difficulties  with  this  approach  for  the  specification  of  a  are  (1) 
the  qualitative  nature  of  the  descriptive  relationship  for  a  (e.g. ,  "normal" 
wet  grass,  crops  with  "sufficient"  moisture),  and  (2)  the  lack  of  time 
dependence  to  account  for  changes  in  a  due  to  specific  precipitation  events 
or  drying-out  periods.  As  part  of  the  current  study,  a  simple  scheme  for 
predicting  time-dependent  surface  moisture  from  routlnely-available 
precipitation  data  and  land  use  characteristics  has  been  developed  and  is 
described  in  Section  2.2.2. 

Figure  8  illustrates  the  relationship  between  the  Bowen  ratio  and  the 

surface  moisture  parameter  at  a  temperature  of  20°C  for  four  values  of  AQ 

2 

(i.e. ,  50,  100,  300,  and  600  W/m  ).  Typical  values  of  the  Bowen  ratio  for 
urban,  suburban,  and  rural  areas  are  1.5,  1.0,  and  0.5,  respectively.  Figure 
7  indicates  that  these  values  of  Bowen  ratio  correspond  to  values  of  a  of 
approximately  0.55,  0.7,  and  0.9,  respectively,  for  typical  mid-day  values  of 
AQ  (i.e.,  300  to  600  W/m^).  Thus,  the  values  relating  a  to  B  with  Eqn.  (14) 
are  consistent  with  typical  empirical  values  of  a  suggested  by  de  Bruin  and 
Holtslag  (1982),  Holtslag  and  van  Ulden  (1983),  and  Hanna  and  Chang  (1991). 

The  typical  value  and  range  of  Bowen  ratio  given  by  Oke  (1978)  and 

2 

corresponding  values  of  a  as  determined  by  Eqn.  (16)  with  AQ  =  300  to  600  W/m" 
are: 

Urban:  B  =  1.5  (range  0.5  to  over  4.0) 

a  =  0. 55  (range  0.2  to  0.9) 
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Bowen  Ratio 


Suburban: 


B  =  1.0 
a  =  0.7 

Rural:  8  =»  0.  5 

a  *  0.  9 


(range  0.  25  to  2.5) 
(range  0.  35  to  1.0) 

(range  0. 10  to  1.5) 
(range  0.55  to  1.2) 


2.3.1  Soil  Moisture  -  Numerical  Models 


Many  coupled  soil  moisture-surface  temperature  numerical  models  have 
been  developed  to  simulate  sensible  and  latent  heat  fluxes  from  bare  soil  or 
vegetated  surfaces.  Although  the  complexity,  high  computational  cost,  and 
extensive  data  input  requirements  of  these  multi-layer  prognostic  models  make 
them  unsuitable  for  use  In  the  current  study,  the  results  from  the  complex 
modeling  studies  are  useful  in  understanding  the  behavior  of  soil  moisture 
variability  and  can  be  helpful  in  developing  simple  parame ter izat ions  of  the 
behavior  of  soil  moisture. 


Deardorff  (1977)  developed  a  model  for  evaporation  over  bare  soil 

surfaces  which  included  prognostic  equations  for  soil  moisture  which 

distinguished  between  soil  surface  moisture  and  bulk  soil  moisture  content. 

The  evaporation  rate  is  parameterized  in  terms  of  the  surface  volumetric 

water  content,  w  ,  (reflecting  conditions  in  the  top  few  millimeters  of  the 
S 

soil)  rather  than  the  bulk  volumetric  water  content,  w^,  (representing 
average  volumetric  soil  moisture  through  a  deeper  layer  from  the  surface  down 
to  about  50  cm). 
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where  wgat  is  the  soil-surface  moisture  content  above  which  evaporation 

occurs  at  the  potential  rate  (w  .  -  0.75  w  ),  and 

sat  max 

Epot  is  the  potential  evaporation  rate. 

Thus,  the  addition  of  moisture  to  the  system  when  the  surface  moisture 

content  is  at  or  above  0.75  w  does  not  increase  surface  evaporation  rate. 

max 

in  Deardorff ’ s  model,  the  surface  moisture  varies  in  response  to  the 
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difference  in  surface  evaporation  and  precipitation  and  includes  a  term 
tending  to  exponentially  restore  surface  moisture  to  the  bulk  value  (w,q)  in 
the  absence  of  significant  evaporation  or  precipitation.  However,  the  model 
does  not  treat  drainage  from  the  thick  (bulk)  soil  layer  or  runoff  losses. 

In  addition,  because  the  bulk  soil  layer  rather  than  the  surface  layer 
supplies  the  moisture  for  transpiration  from  vegetation,  this 
parameterization  is  most  appropriate  for  bare  soil  surfaces.  However, 
Deardorff  (1978)  presents  a  generalized  model  which  produces  separate 
estimates  for  transpiration  rates  in  vegetated  areas  and  surface  evaporation 
over  areas  without  vegetation. 

However,  one  of  the  significant  results  of  the  Deardorff  (1977)  model  is 

the  important  effect  of  rapid  soil  surface  drying  on  reducing  the  rate 

of  evaporation.  Rapid  surface  drying  is  expected  to  be  important  in 

urban  environments,  in  which  a  large  percentage  of  the  surface  is  covered 

with  impervious  surfaces,  and  drainage  systems  exist  to  handle  runoff. 

Figure  9  contains  predicted  values  of  w^  and  w^  using  Deardorff’ s  model 

and  observed  values  of  w  for  a  field  of  bare  soil.  The  observed  data  is 

8 

from  Jackson  (1973)  for  a  field  near  Phoenix,  Arizona  during  March  1971  for  a 
seven  day  period  following  irrigation  of  the  soil.  Deardorff  (1977)  notes 
that  the  use  of  the  bulk  moisture  content  would  have  overpredicted 
evaporation  rates  by  a  factor  between  1.5  and  5.6  during  days  6  and  7. 

Acs  et  al.  (1991)  developed  a  three-layer  soil  moisture  and  soil  surface 
temperature  model.  The  volumetric  soil  moisture  content  is  based  on  the 
method  of  Sellers  et  al.  (1986). 
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Figure  9. 


Predicted  values  of  surface  volumetric  water  content,  w  (dashed 

S 

line),  and  bulk  volumetric  water  content,  w^  (dotted  line).  Also 
shown  are  observed  values  of  w  (solid  line),  predicted 

s 

evaporation  (mm)  and  observed  evaporation  (mm)  in  parentheses. 
(From  Deardorff,  1977;  observed  data  from  Jackson,  1973). 
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where  w^  is  the  volumetric  soil  moisture  of  the  ith  layer, 

P  is  the  infiltration  due  to  precipitation, 

E  is  the  soil  surface  evaporation  rate, 
is  the  thickness  of  the  ith  soil  layer, 

Pw  is  the  density  of  water, 

^+1  is  the  flow  between  soil  layer  1  and  1+1,  and 
is  the  drainage  of  water  from  the  bottom  of  soil  layer  i. 

The  layers  in  Acs  et  al.  (1991)  model  are  0-10  cm,  10-30  cm,  and 
30-60  cm  in  depth.  Detailed  information  on  the  hydraulic  properties  of  the 
soil,  such  as  soil  pore  space  and  hydraulic  conductivity  at  saturation,  are 
used  to  compute  the  water  flow  between  layers.  The  evaporation  rate  is 
assumed  to  be  equal  to  the  potential  rate  when  w^  exceeded  0.  75  wsat  where 
wsat  is  the  volumetric  soil  moisture  at  saturation,  and  a  linear  function  of 
w1  when  Wj  is  below  0.75  «sat-  Figure  10  Illustrates  the  strong  diurnal 
variation  of  soil  moisture  in  the  surface  layer.  Pronounced  daytime  drying 
is  predicted  along  with  increases  in  soil  moisture  during  nighttime  periods , 

although  the  model  tends  to  overpredicted  drying  during  the  5-day  periods. 


Pan  and  Mahrt  (1987)  investigated  the  interactions  among  soil  moisture, 
energy  fluxes,  and  boundary  layer  development.  They  used  a  detailed  34-layer 
atmospheric  boundary  layer  model  (extending  up  to  4  km  in  height  )  coupled 
with  a  two-layer  soil  model  (5  cm  deep  top  layer,  95  cm  deep  bottom  layer)  in 
the  study.  They  describe  three  stages  of  soil  drying  with  different  effects 
on  energy  fluxes  and  boundary  layer  development.  These  stages  are 
illustrated  in  Figure  11.  In  the  first  stage,  surface  evaporation  is  at  the 
potential  rate  and  is  determined  primarily  by  wind  speed,  relative  humidity, 
and  Incoming  solar  radiation.  Surface  evaporation  is  strong,  which  limits 
surface  heating  and  therefore  boundary  layer  development.  The  potential  and 
actual  evaporation  rates  decrease  moderately  with  time  due  to  a  gradually 
moistening  of  the  boundary  layer. 
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toil  motstur?  content  ( TO* Jn**) 


Figure  10,  Predicted  soil  moisture  in  three  layers  [0-10  cm  (solid  line), 

10-30  cm  (fine  dashed  line),  and  30-60  cm  (coarse  dashed  line)]. 
Periods  with  precipitation  are  indicated  with  arrows.  Data  for  a 
surface  of  bare  soil  during  15-20  April,  1988.  (From  Acs  et  al. , 
1991). 
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Figure  11.  Predicted  evolution  of  (a)  noontime  surface  latent  heat  flux 
2 

(W/ra  ),  (b)  noontime  boundary  layer  height,  and  (c)  noontime 
potential  evaporation  for  two  surface  types  (clay  and  sand) 
during  periods  of  moderate  geostrophic  winds  (5  m/s)  and  strong 
geostrophic  winds  (10  m/s)  on  the  summer  solstice.  (From  Pan  and 
Mahrt,  1991). 
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In  the  second  stage,  the  evaporation  rate  decreases  rapidly  below 
potential  evaporation  values  due  to  drying  of  the  surface  layer  (Figure 
11(a)).  Actual  evaporation  is  limited  by  the  rate  of  soil  moisture  transport 
to  the  surface,  which  is  determined  by  soil  properties  and  the  moisture 
gradient  within  the  soil.  With  a  lower  surface  evaporation  rate,  surface 
temperature,  sensible  heat  flux  and  boundary  layer  height  all  increase 
significantly  (Figure  11(b)).  In  the  second  stage,  the  entrainment  of  drier 
air  from  above  the  boundary  layer  can  be  significant,  ev  n  exceeding  the 
surface  evaporation,  and  leading  to  an  overall  drying  of  the  boundary  layer. 
Even  though  potential  evaporation  increases  due  to  higher  temperatures  and 
boundary  layer  drying  (Figure  11(c)),  the  actual  evaporation  rate  decreases 
due  to  limited  surface  moisture  availability. 

The  final  stage  of  drying  occurs  near  steady-state  conditions  are 
reached.  Eventually  the  mid-day  surface  evaporation  rat*  decreases  to  a 
small  fraction  of  the  potential  rate  (i.e.  ,  <  1054).  The  soil  moisture  is 
significantly  cepleted  and  changes  further  only  slowly.  The  boundary  layer 
development  is  deep,  and  is  characterized  by  warm  and  dry  condition-. 

Pan  and  Mahrt  (1987)  conducted  sensitivity  tests  of  the  modeling 
results  on  factors  such  as  vegetation  and  solar  radiation.  Vegetation  tends 
to  reduce  the  differences  between  ‘he  three  stages.  Because  vegetation  draws 
moisture  from  a  deeper  layer  (i.e.,  the  root  zone),  soil  drying  occurs  more 
slowly  and  the  onset  and  magnitude  of  the  stage  2  effect  is  reduced  (see 
Figure  12(a)).  Reduced  solar  radiation  lowers  the  potential  evaporation  rate 
which  also  slows  the  drying  process.  At  wintertime  levels  of  solar 
radiation,  the  vertical  transport  of  coil  moisture  from  deeper  soil  layers 
during  nighttime  periods  can  restore  or  nearly  restore  the  surface  moisture 
level  sufficiently  to  maintain  daytime  potential  evaporation  rates  (see 
Figure  12(b)). 

Noilhan  and  Planton  (1989)  used  a  force  restore  method  based  on 
Deardorff  (1977,  1978)  to  parameterize  surface  moisture  and  its  effect  on 
latent  and  sensible  heat  fluxes.  The  volumetric  soil  moisture  content  is 
described  by: 
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(a) 


Day 


(b) 


Figure  12.  Results  of  sensitivity  tests  showing  (a)  noontime  surface  latent 

2 

heat  flux  (W/m  )  for  a  vegetated  surface  on  the  summer  nolstice, 

2 

and  (b)  noontime  surface  latent  heat  f^ux  (W/m  )  for  a  bare 
surface  on  the  winter  solstice.  (From  Pan  and  Mahrt,  1991). 
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where  w  is  the  volumetric  soil  moisture  content  in  a  thin  surface  soil  layer 

s 

(few  millimeters  in  depth), 

*2  is  the  average  volumetric  soil  moisture  content  through  a  deeper 
layer  from  the  surface  down  to  a  depth  d^  (~  50  cm), 
w  is  the  equilibrium  surface  volumetric  soil  moisture  content  when 
gravity  balances  the  capillarity  forces, 

Et  is  the  transpiration  rate, 
t  is  a  time  constant, 

dj  is  a  normalization  constant  (~  10  cm),  and, 

Cj,  are  empirical  coefficients. 


The  other  variables  have  been  defined  previously.  In  their  model, 

the  atmospheric  conditions  forcing  the  system  (i.e. ,  solar  radiation  wind 

speed,  specific  humidity,  temperature)  are  specified  rather  than  computed. 

The  most  important  improvements  of  Deardorff’s  model  made  by  Noilhan  and 

Planton  (1989)  are  (1)  a  correction  to  Deardorff’s  restore  term  in  Eqn.  (22) 

to  account  for  gravitational  effects  (i.e.,  the  use  of  w  instead  of  w_  in 

geq  2 

the  last  term  of  the  equation)  and  (2)  calibration  of  the  coefficients  in  the 
equation  based  for  different  types  and  soil  wetness  values. 


2.3.2  Soil  Moisture  -  Capacitance  Model 


A  simple  capacitance  model  has  been  developed  to  allow  the  soil  moisture 
parameter,  a,  to  be  estimated  operationally  from  routinely-available 
precipitation  data  and  land  use  information.  The  model  is  intended  to 


37 


reproduce  some  of  the  most  significant  features  of  the  behavior  of  the 
comprehensive  soil  moisture-surface  evaporation  models  discussed  in  the 
previous  section,  but  with  significantly  more  modest  input  and  computational 
requirements. 


The  model  consists  of  two  capacitors:  a  slow-response  capacitor  to 
simulate  bulk  moisture  content  through  a  deep  layer  and  a  fast-response 
capacitor  to  respond  to  moisture  variations  in  a  thin  surface  layer. 
Precipitation  acts  to  "charge"  the  capacitor  and  evaporation  “discharges"  the 
system.  Feedback  is  allowed  between  the  bulk  and  surface  moisture  parameters 
to  simulate  a  build-up  ("re-charge")  of  surface  moisture  due  to  transport 
from  deeper  layers  during  periods  of  low  evaporation  (e.g. ,  see  Figure  9). 
Different  values  of  the  capacitor’ s  time  constant  are  used  to  parameterize 
the  variable  responses  to  moisture  input  of  different  surface  (land  use) 
types  (e.g.,  urban,  suburban,  agricultural,  vegetated  surfaces).  In 
converting  the  "charge"  (q)  to  a  soil  moisture  parameter  (a),  a  cap  is 
applied  to  a  which  results  in  a  nearly  constant  evaporation  rate  in  the 
initial  (stage  1)  phase  of  drying  from  a  moist  surface,  followed  by  a  period 
of  rapid  decrease  in  evaporation  (l.e. ,  the  stage  2  effect),  and  finally, 
very  low,  near  steady-state  evaporation  rates  characteristic  of  stage  3  fsee 
Figure  11 (a) 1 


During  drying-out  periods. 


q(t+At)  *  q(min)  +  Aq  exp  -< 


Aq  *  q(max)  -  q(rain) 


the  state  of  charge  of  a  capacitor 


-(t  +  At) 


t  at;  ^ 

~  f 


is: 


(24) 

(25) 


where  q  is  the  “charge"  of  the  system  at  time  (t+At), 

q(min),  q(max)  are  minimum  and  maximum  values  of  q  for  a 
particular  land  use  type, 

r  is  a  time  constant  at  time  (t+At)  which  depends  on  land  use  type 
and  time  of  day,  and, 

t  is  a  virtual  time  computed  from  q  at  the  previous  time  step 
(i.  e. ,  q ( t ) ) . 
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The  time  constant  and  minimum/maximum  values  of  q  must  be  determined 
empirically.  Additional  work,  beyond  the  scope  of  the  current  study,  is 
required  to  develop  tables  of  '‘typical'*  values  of  these  parameters  as  a 
function  of  routinely-available  data,  such  as  land  use.  However,  because  the 
system  is  likely  to  be  heavily  driven  by  the  measured  precipitation  history, 
it  is  expected  that  the  results  will  not  be  overly  sensitive  to  errors  in  the 
specification  of  the  empirical  parameters  r,  q(min),  and  q(max). 

It  is  necessary  to  compute  a  virtual  time,  t  ,  because  r  may  vary 
from  one  time  step  to  the  next.  Therefore,  t  is  computed  to  reflect  the 
previous  state  of  the  system  in  terms  of  the  new  r  at  time  t+At. 


-r  In 


q(t)  -  q(min) 
Aq 


} 


(26) 


When  the  system  is  being  recharged  (i.e. ,  during  a  precipitation  event), 
the  state  of  the  system  is: 


q(t+At) 


q(min)  +  Aq 


{  1  -  exp  [ 


-  (t 


At) 


]} 


(27) 


where  the  virtual  time  is: 


In 


q(t)  -  q(min) 

_ 


} 


(28) 


The  time  constant  during  precipitation  periods  is  computed  from  the 
precipitation  amount: 


t  -  C/R 


(29) 


where  £  is  an  empirical  parameter  (hour/mm)  which  depends  on  the 
surface  type,  and, 

R  is  the  precipitation  amount  (mm). 


The  surface  soil  moisture  parameter,  a,  is  estimated  from  the  capacitor 
charge,  q,  as: 
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a  =  •< 


q 


q  <  1.0 


(30) 


(1.0  q  >  1.0 

Thus,  the  addition  of  moisture  to  the  system  when  it  Is  at  or  beyond  fully 
moist  conditions  does  not  change  the  soil  moisture  parameter  or  the  surface 
evaporation  rate.  Also,  the  surface  evaporation  rate  will  remain  nearly 
constant  during  the  initial  phase  (i.e. ,  stage  1)  of  a  drying-out  period 
until  q  drops  below  a  value  of  one.  However,  q  (and  a)  will  vary  rapidly 
with  continued  drying  once  q  decreases  to  be  less  than  one  (i.e.,  stage  2). 
Further  drying  will  result  in  q  approaching  a  constant  value  of  q(min)  (i.e., 
stage  3). 

In  the  capacitance  model  used  in  this  study,  two  values  of  q  are 
tracked:  q^.  is  the  value  for  the  fast-response  capacitor  representing 
moisture  in  a  thin  surface  layer  and  qg  is  the  value  for  the  slow-response 
capacitor  (tracking  bulk  moisture  content).  Each  capacitor  has  its  own  value 
of  £  and  t  (i.e.,  <;  and  r_,  t  )  reflecting  the  different  rates  of 

response  to  precipitation  Input  and  evaporation  output.  Additional 
capacitors  with  slower  response  time  simulating  deeper  (or  thicker)  layers 
could  be  added.  However,  the  two-capacitor  model  is  used  to  demonstrate  the 
technique  and  appears  to  be  adequate.  Although  t  Is  a  function  of  the 
evaporation  rate,  no  attempt  has  been  made  to  compute  hourly  values  for  the 
current  set  of  tests,  although  this  could  be  a  future  enhancement.  However, 
different  values  of  t  have  been  used  for  daytime  and  nighttime  periods  to 
reflect  the  large  differences  in  evaporation  rates  between  these  periods. 

During  nighttime  periods  when  evaporation  rates  are  low,  the  surface 
moisture  content  tends  to  be  restored  toward  the  bulk  value  due  to  moisture 
transport  from  deeper  layers  to  the  surface.  If  the  nighttime  surface 
moisture  is  greater  than  the  >ulk  value,  the  surface  moisture  will  tend  to 
slowly  decrease  and  approach  the  bulk  value  from  above'.  These  features  of 
the  soil  moisture  behavior  are  reproduced  by  setting  q(min)  -  qg  (if  q^.  >  qg) 
and  q(max)  3  qg  (if  q^.  <  q  ). 

The  daytime  drying  and  nighttime  restoring  of  the  surface  moisture 
predicted  oy  the  capacitance  model  is  shown  in  Figure  13  for  a  seven-day 
drying  period.  The  assumed  values  for  the  model  parameters  are 
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Surface  Moisture  Parameter 


Figure  13.  Illustration  of  daytime  drying  and  nighttime  moistening  of 

surface  soil  as  predicted  by  the  capacitance  model  for  a  seven-day 

drying  out  period.  Fast-response  (surface  layer)  value  (q^J  and 

slow-response  (bulk  layer)  value  (qg)  of  soil  moisture  are 

plotted.  Model  parameters  are  xf(day)  =  xf (night)  =  15  hr, 

r  (day)  =  80  hr.  and  x  (night)  =  99999  hr. 
s  s 
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T^(day)  =  Tonight)  *  15  hr,  Tg(day)  =  80  hr,  and  Tonight)  =  99999  hr. 

The  predicted  surface  moisture  cycle  shown  in  Figure  13  can  be  qualitatively 
compared  to  the  observed  behavior  shown  in  Deardorff  (1977)  for  Days  3  and 
5-7  (Figure  9). 

The  capacitance  model  has  been  used  estimate  soil  moisture  for  the 
base-case  conditions  described  in  the  21-day  simulation  of  a  drying-out 
period  by  Pan  and  Mahrt  (1991),  The  soil  moisture  is  predicted  with  the 
capacitance  model  for  two  surface  types  (sand  and  clay)  using  the  following 
assumed  model  parameters: 

clay:  Tg(day)  *  65  hrs  sand:  rg(day)  *  20  hrs 

rg (night)  »  99999  hrs  Tg (night)  *  99999  hrs 

Assuming  noontime  conditions  of  AQ  *  600  W/m  ,  T  *  25°C,  and  £'*20  W/m  , 

2 

and  using  Eqn.  (14),  the  latent  heat  flux  (in  W/m  )  can  be  expressed  as: 

Q  (noon)  ■  464  a  (31) 

e 

where,  a  is  determined  using  Eqn.  (30)  with  q=qg.  Figure  14  shows  the 
predicted  noontime  latent  heat  flux  using  soil  moisture  predicted  by  the 
capacitance  model  along  with  the  predictions  from  the  numerical  model  of  Pam 
and  Mahrt  (1991).  The  capacitance  model  successfully  reproduces  the  main 
features  of  the  more  complex  model  for  both  surface  types.  Each  of  the  three 
stages  of  evaporation  discussed  above  can  be  seen  in  the  predictions  of 
both  models. 

The  relationship  between  a  and  q  depends  on  the  type  of  surface  being 
modeled.  For  example,  as  found  by  Deardorff  (1977)  and  (1978),  evaporation 
from  a  bare  soil  surface  is  determined  by  surface  moisture  content.  However, 
over  a  vegetated  surface,  moisture  from  the  deeper  (bulk)  layer  controls 
evapotranspiration.  Therefore,  in  that  case,  a  is  most  appropriately  related 
to  qg.  In  our  testing  of  the  model  for  an  urban  environment,  a  short  time 
scale  was  assigned  to  t^.  to  allow  q^.  to  represent  moisture  on  the  surface  of 
the  urban  surfaces  (e.g.,  wet  pavement,  wet  ground  surfaces).  The  time  scale 
for  Tg  was  selected  to  reflect  longer  term  average  urban  moisture  conditions. 
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Latent  Heat  Flux  (W/m2) 


700 


2 

Figure  14.  Predicted  noont’me  surface  latent  heat  flux  (W/m  )  based  on  soil 
moisture  predictions  of  capacitance  model  for  clay  soil  (•)  and 
sandy  soil  (o)  for  a  21 -day  drying  out  period.  Also  shown  are 
the  numerical  model  predictions  of  Pan  and  Mahrt  (1991): 
clay  (short  dash),  clay,  windy  conditions  (long  dash), 
sand  (intermediate  dash),  sand,  wind  conditions  (solid). 
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The  moisture  parameter  a  was  determined  by  applying  Eqn.  (30)  to  the  larger 
of  q^.  or  qg.  Thus,  a  Increases  rapidly  with  the  occurrence  of  precipitation 
as  the  urban  surfaces  become  wet,  but  relatively  quickly  returns  to  its 
longer  term  value  determined  by  qg  after  the  surface  moisture  evaporates. 

As  an  illustration  of  the  capacitance  model  technique,  predicted  values 
of  q  for  a  three  month  period  are  plotted  in  Figure  15(a)  and  16(a)  for  an 
urban  area  and  vegetated  rural  area,  respectively.  Daily  precipitation 
values  for  Boston,  Massachusetts  during  the  period  from  May  through  July, 

1989  were  used  in  both  simulations  (see  Figure  15(b)).  Monthly  total 
precipitation  for  May  and  June  were  near  climatological  mean  values  (89.9  mm 
and  72.1  mm,  respectively,  representing  deviations  from  normal  of  +0.5  mm 
and  -2.0  mm).  Precipitation  for  July  (129.3  mm)  was  above  normal  by  61.2  mm. 

The  two-capacitor  moisture  model  was  run  for  92  days  starting  on  May  1 

with  an  initial  value  of  q  of  0.5  for  both  areas.  The  main  difference  in  the 

urban/ rural  inputs  assumed  were  the  values  of  £g  (96  hr/mm  (urban)  vs.  48  hr/mm 

(rural))  and  r  (24  hr  (urban)  vs.  96  hour  (rural)).  The  larger  value  of  £ 
s  s 

for  the  urban  area  reflects  the  diversion  of  a  fraction  of  the  precipitation 
from  soil  storage  in  the  urban  environment  by  impervious  surfaces,  urban 
drainage  systems  and  runoff.  Thus,  a  given  aunount  of  precipitation  in  an 
urban  area  will  produce  a  smaller  change  in  the  bulk  soil  moisture  content 
than  in  a  rural  area  where  more  of  the  precipitation  is  absorbed  into  the 
ground.  Similarly,  values  of  Tg  were  selected  to  allow  a  more  rapid  drying 
of  the  urban  area.  Unlike  the  rural  environment,  where  the  precipitation 
makes  its  way  into  surface  and  deeper  soil  layers  and  therefore  is  available 
to  be  evaporated  over  the  next  few  days,  much  of  the  precipitation  input  into 
the  urban  system  is  carried  away,  resulting  in  a  more  rapid  return  to  normal 
dry  conditions.  Similar  values  of  (12  hrs  (urban),  15  hrs  (rural))  and  ^ 
(24  hr/mm  for  both)  were  used  in  the  simulations. 

The  predictions  of  the  capacitance  model  in  Figure  15(a)  are  generally 
consistent  with  observations  of  the  Bowen  ratio  (related  to  a  by  Eqn.  (16)) 
as  described  in  Section  2.2.  Typical  observed  urban  values  of  the  Bowen 
ratio  are  ~  1.5,  which  corresponds  to  a  mid-day  a  of  -  0.5.  (According  to 
Eqn.  (30),  a=q  when  q  <  1.0,  and  a  =  1.0  when  q  >  1.0).  The  predicted  daily 
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averaged  urban  value  of  a  for  most  af  the  days  without  precipitation  are  in 
the  range  0.4  to  0.6  (i.e. ,  on  40  of  55  dry  days)  with  am  average  dry  day 
value  of  0.5.  During  precipitation  events,  q  (and  a)  increase  rapidly  as  a 
result  of  precipitation  wetting  of  the  urban  surface.  However,  a  also  falls 
quickly  after  the  precipitation  event,  reflecting  rapid  drying  of  the  urban 
environment.  In  Figure  16(b),  the  soil  moisture  parameter  in  the  rural  area 
stays  within  the  range  0.8  to  1.0  for  most  of  the  period.  These  values  of  <x 
are  consistent  with  observed  rural  Bowen  ratios  ~  0. 5.  During  a  six-day 
period  with  no  precipitation  in  late  June  and  early  July,  drying  results  in  a 
drop  in  a  to  0.68.  Also  note  that  the  low  initial  value  of  a.  *  0. 5  is 
quickly  adjusted  by  the  model  to  more  typical  rural  values. 

The  results  discussed  in  this  section  suggest  that  the  capacitance  model 
may  be  a  useful  tool  in  estimating  the  time  variability  of  the  soil  moisture 
parameter,  which  is  one  of  the  key  inputs  to  the  energy  balance  model.  The 
capacitance  model  produces  reasonable  values  of  a  for  both  urban  and  rural 
environments.  It  requires  as  input  only  routinely-available  precipitation 
data  and  empirical  time  scale  parameters  (i.e.,  t  and  £)  which  can  be  related 
to  land  use.  The  test  case  simulations  produce  estimates  of  soil  moisture 
which  are  consistent  with  observations  of  typical  Bowen  ratios  in  urban  and 
rural  areas. 

2.  4  Parameterization  of  Boundary  Layer  Parameters 

A  number  of  meteorological  preprocessors  have  been  developed  which  employ 
an  energy  balance  method  based  on  Holtslag  and  van  Ulden  (1983)  to  estimate 
the  sensible  heat  flux  from  routinely  available  meteorological  data  and  then 
compute  boundary  layer  parameters  such  as  the  surface  friction  velocity  and 
Monin-Obukhov  Length  with  empirical  or  iterative  schemes.  Such  preprocessors 
include  METPRO  (Hanna  et  al. ,  1986),  SIGMET  (Scire  and  Wojichowski,  1987), 
SIGPRO  (Hanna  and  Chang,  1990  and  1992)  and  CALMET  (Scire  et  al. ,  1990a),  as 
well  as  several  others.  Of  these,  all  except  METPRO  contain  features  which 
allow  them  to  be  used  In  urban  areas.  However,  because  SIGPRO  has  several 
enhanced  features  and  has  the  advantage  of  being  recently  evaluated  using 
data  from  urban  field  experiments  in  St.  Louis  and  Indianapolis,  it  is 
recommended  as  the  starting  point  for  the  Phase  II  study. 
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SIGPRO  uses  the  sensible  heat  flux,  Q^,  produced  from  the  energy  balance 
model  (see  Section  2.2)  along  with  observations  of  the  wind  speed,  u,  and 
estimates  of  the  surface  roughness  length,  zqI  and  displacement  height,  d,  to 
compute  the  following  boundary  layer  parameters: 

u„  -  surface  friction  velocity  [ (-u7wT)1/2] , 

L  -  Monin-Obukhov  length, 
h  -  mixing  height,  and, 
w„  -  convective  velocity  scale 


Because  SIGMET  allows  the  input  of  hourly  values  of  the  surface  moisture 
parameter,  a,  it  cam  be  easily  coupled  to  the  capacitance  soil  moisture  model 
described  in  Section  2. 3. 2. 


2.4.1  Friction  Velocity  and  Monin-Obukhov  Length 

The  friction  velocity  is  expressed  in  terms  of  the  following  wind 
profile  formula: 


u 


-  (z/L) 

in 


} 


(32) 


where  u  is  the  wind  speed  at  height  z, 

k  is  the  von  Karman  constant  (~0.  4), 

zq  is  the  surface  roughness  length,  and, 

is  an  universal  stability  correction  factor  which  is  a 
m 

function  of  z/L  (e.g. ,  Dyer  and  Hicks,  1970). 


Hanna  and  Chang  (1992)  recommend  the  following  rules  of  thumb  for 

d  and  z  :  d  ~  0. 5h  and  z  -  0. lh  ,  where  h  is  the  height  of  roughness 
o  r  o  r  r  °  ° 

elements  (Stull,  1988).  The  Monin-Obukhov  length  is  defined  as: 


L 


3  T 

-  u.  T  p  c 
*  P 

k  gQh 


(33) 
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where  T  is  the  air  temperature, 
p  is  the  air  density, 

c  is  the  specific  heat  at  constant  pressure, 
P 

g  is  the  acceleration  due  to  gravity. 


and, 


If  the  buoyancy  effects  of  water  vapor  are  to  be  included,  in 
Eqn.  (33)  is  replaced  with  the  virtual  heat  flux,  Q^v,  and  T  is  replaced  with 
the  virtual  temperature,  T  ,  as  described  in  Section  2. 4. 2. 

NEUTRAL  CONDITIONS 


Under  neutral  conditions,  i <j)  is  zero,  and  u.  can  be  computed  directly 
from  Eqn.  (32),  i.e. , 


_ ku _ 

u*n  ~  ln[(z  -d)/z  ] 
o 


(34) 


Once  u„  is  known  and  is  determined  from  the  energy  balance  equations, 
L  is  computed  directly  from  its  definition  (Eqn.  (33)). 

UNSTABLE  CONDITIONS 

When  is  not  zero,  solution  of  Eqns.  (32)  and  (33)  requires  the 
application  of  a  computationally  demanding  iterative  technique.  An 
alternative  approach  used  in  SIGPRO  and  SIGMET  for  unstable  conditions  is  an 
empirical  relationship  developed  by  Wang  and  Chen  (1980)  which  can  be  solved 
analytically. 

u*  *  lnllf  -  d)/20]  ['  *  dlln(1  *  W]  1351 

where 

d.  =»  0.128  +  0. 0051n(z  /z)  if  z  /z  <  0.01  (36) 

l  o  o  - 

if  z  /z  >  0.U1 
o 


=  0. 107 


(37) 


d 
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2 
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1.95  +  32.6(2  /z) 
o 


0.  45 


%  kgz 


pc  Tu. 

K  p  *n 


(38) 

(39) 


The  term  d^ln(l  +  d^d^)  represents  the  correction  due  to  instability  and  u„n 
is  defined  in  Eqn  (34). 


Hanna  and  Chang  (1990,  1992)  tested  the  analytical  formula  against  values 

produced  by  the  iterative  solution  of  u.  and  L.  They  found  that  the  Wang  and 

Chen  (1980)  expression  produced  values  within  1054  of  the  results  determined  by 

the  iterative  solution  for  z  =*  10  m,  d  *  0,  z  aim,  and  a  large  value  of  Q. 

2  oh 

(400  W/m  ).  Better  agreement  was  found  for  smaller  roughness  elements  and 

smaller  sensible  heat  fluxes.  In  addition,  the  analytical  solution  was 

computationally  significantly  faster.  Therefore,  Hanna  and  Chang  (1990) 

recommend  the  use  of  the  Uang  and  Chen  (1980)  approach. 

Again,  once  u„  and  are  known,  L  is  computed  directly  from  Eqn.  (33). 

STABLE  CONDITIONS 


The  Weil  and  Brower  (1983)  method  for  estimating  u»  is  applied  in  SIGPR0 
during  stable  conditions.  A  first  estimate  of  the  scaling  temperature,  9,, 
is  calculated  using  Holtslag  and  Van  Olden’ s  (1983)  equation: 


9B1  *  0.09(1  -  0.5N2)  (40 

where  N  is  the  total  cloud  cover  in  fractions*  and  0,  has  units  of  K. 
Another  estimate  of  0,  is  made  from  the  profile  equation  for  temperature: 


0 


•2 


TC  .  u 
dn 

18. Szg 


(41) 


where  the  neutral  drag  coefficient  C,  is  defined  as  k/£n[(z  -  d)/z  ]. 

dn  o 


Then,  9.  is  set  equal  to  the  smaller  of  9,^  and  9m^. 
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The  sensible  heat  flux,  Q.,  is  defined  during  stable  conditions  as: 

n 

\  =  -  pcpu„9..  (42) 

For  large  values  of  u  (or  u„),  9.^  (which  depends  only  on  cloud  cover)  is 
smaller  than  9, 2>  but  am  additional  check  on  the  product  u.9.  must  be  made, 
since  does  not  keep  increasing  indefinitely  with  higher  wind  speeds.  In 
SIGPRO,  the  value  of  9,  is  not  allowed  to  exceed  Q.05/u,,  where  the 
numerator  has  units  of  ta/s  and  denominator  has  units  of  m/sec.  This  limit 
is  estimated  from  observations  of  neat  fluxes  during  high-wind,  stable 
conditions. 


The  friction  velocity,  u„,  can  be  calculated  from: 


u„  - 


C  .  u 
dn 


1  ♦ 


1  - 


2u  \  ‘ 

n7s — ) 

dn 


2\  1/2 


(43) 


where  uq  =  ^4. 7zg9,/T| 


1/2 


Because  9.  is  set  equal  to  the  smaller  of  9#1  and  9»2,  the  following 
condition  is  always  met: 


2u 


„  1/2  - 
C ,  u 
dn 


<  1 


(44) 


During  stable  conditions,  SIGPRO  enforces  a  lower  limit  on  L  in 
recognition  of  the  fact  that  the  atmosphere  is  less  stable  over  urban  areas 
than  over  rural  surfaces.  The  minimum  values  of  L  used  as  default  values  in 
SIGPRO  are  based  on  the  observation  that  the  mechanical  mixed  layer  is  about 
2  to  3  times  the  typical  building  height  (Uno  et  al.  ,  1988)  and  that  L 
represents  the  height  of  the  mechanically  mixed  layer.  Hanna  and  Chang  (1990) 
suggest  the  values  in  Table  4  for  use  for  the  various  land  use  categories 
defined  in  the  Auer  (1978)  scheme. 


51 


Table  4 


Minimum  Vaiues  of  Mon in-Obukhov  Length 
During  Stable  Conditions 
for  Various  Land  Use  Types 
(From  Hanna  and  Chang,  1990) 


Auer  (1978) 

Class 

Description 

Minimum  L 

Cl 

Commercial 

>  40  story  buildings 

150  m 

10-40  story  buildings 

100  m 

<  10  story  buildings 

50  a 

11, 

12 

Industrial 

50  m 

R3 

Compact  Residential 

50  m 

Rl, 

R2 

Residential 

25  a 

A 

Agricultural 

2  m 
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2.4.2  Inclusion  of  Buoyancy  Effects  of  Water  Vapor 


During  most  conditions  over  land  surfaces,  evaporation  rates  are 
suff iciently  small  that  the  buoyancy  effects  of  water  vapor  are  negligible  In 
the  expressions  for  the  boundary  layer  parameters.  However,  evaporation  rates 
can  sometimes  be  high  enough  (e.g. .  after  a  summertime  precipitation  event 
with  high  net  radiation  conditions)  that  water  vapor  can  affect  the  density 
stratification  near  the  surface.  The  effects  of  water  vapor  buoyancy  can  be 
accounted  for  by  using  a  virtual  heat  flux,  Q^v  Instead  of  and  a  virtual 
temperature  (Tv)  instead  of  T  in  the  relevant  equations  for  L  (Eqn.  (33))  and 
w,  (Eqn.  (50)).  Arya  (1988)  suggests  the  following  formula  to  estimate  Q^: 

+  °-61  cp  9Qe/Le  (45) 

where  9  is  the  mean  potential  temperature, 

3 

c  is  the  specific  heat  at  constant  pressure  (10  J/Ocg  K)), 

P  g 

Le  is  the  latent  heat  of  evaporation  (2.5x10  J/kg),  and, 

Q  is  the  latent  heat  flux, 
e 

The  virtual  temperature  is  the  temperature  of  dry  air  having  the  same 
density  as  moist  air  at  the  same  pressure.  It  is  approximated  as: 

Ty  ~  T  (1  +  0.  6  q^)  (46) 

where  q  is  the  mixing  ratio  (q  =  RH  q  /100. ), 
m  ms 

RH  is  the  relative  humidity  of  the  air,  and 

qs  is  the  saturation  mixing  ratio  (temperature  and  pressure  dependent). 

If  the  data  to  compute  Ty  Is  not  available,  its  effect  can  be 

neglected,  since  the  most  important  effect  of  water  vapor  buoyancy  on  L 

is  accounted  for  in  the  Q^v  terra.  Arya  (1988)  notes  that  Ty  usually  does  not 

vary  from  T  by  more  than  a  few  degrees,  so  g/Tv  ~  g/T.  However,  Q^v  can 

sometimes  be  significantly  larger  than  Q^.  For  example,  at  9  =  280  K,  Eqn. 

(45)  yields  Q^v  ~  +  0.07  or  equivalently  -  Q^( 1  +  0.07/B),  where  B 

is  the  Bowen  ratio.  Since  the  Bowen  ratio  can  reach  values  as  low  as  0.1, 

under  these  conditions  Q,  Is  increased  by  7(3%  over  Q.  . 

hv  h 
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2.4.3  Mixing  Height 


The  convective  boundary  layer  (CBL)  is  assumed  to  be  capped  by  a 
relatively  thin  interfacial  layer  separating  it  from  the  stable  air  aloft. 
Stable  air  is  entrained  into  the  interfacial  layer  as  a  result  of  vertical 
mixing  due  to  thermals  penetrating  the  top  of  the  CBL.  Within  the  mixed 
layer  (0.  lh  <  z  <  h),  the  potential  temperature  9  and  wind  speed  u  are 
relatively  uniform,  but  in  the  interfacial  layer,  they  rapidly  adjust  with 
height  to  their  values  in  the  overlying  stable  air. 

Carson  (1973)  has  simplified  the  modeling  of  boundary  layer  evolution  by 
ignoring  radiation,  latent  heat  effects,  and  advectlon  of  energy.  He 
includes  the  effects  of  1)  time-dependent  surface  heating;  2)  capping  layer 
stability;  3)  large-scale  air  subsidence;  and  4)  a  uniform  distribution  of  9 
with  z  within  the  CBL,  a  step  change  A0^  at  the  top  of  the  CBL  (z  =  h),  and  a 
linear  variation  with  z  in  the  overlying  stable  air.  Carson  also  assumes 
that  the  stable  air  is  composed  of  as  many  as  three  vertical  layers,  each 
with  a  different  lapse  rate,  instead  of  a  single  stable  layer. 

Weil  and  Brower  (1983)  modified  Carson’s  model  by: 

•  permitting  the  elevated  stable  layer  to  have  an  arbitrary 
temperature  distribution  with  z  (l.e. ,  an  infinite  number  of 
vertical  layers ) ; 

•  accounting  for  surface  stress- induced  (mechanical)  mixing,  which 
can  be  important  at  night  and  in  the  early  morning  hours  when  the 
heat  flux  is  low  or  zero;  and 

•  neglecting  subsidence. 

Convective  and  mechanical  mixing  are  assumed  to  be  independent  of  one 
another,  so  that  when  one  of  these  mixing  modes  is  operative,  the  other  is 
not.  Figure  17  shows  an  example  of  an  assumed  potential  temperature 
distribution,  where  the  solid  curve  is  the  initial  temperature  profile, 

9^ ( z ) ,  and  the  dashed  curve  is  the  idealized  profile  at  a  later  time  t.  The 
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Figure  17.  Mixing  height  computation  using  the  modified  Carson  method. 
(From  Hanna  and  Chang,  1990). 
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"overshoot"  <r  is  a  measure  of  the  degree  of  entrainment  or  interfacial  mixing 
and  is  a  function  of  time,  as  are  the  mixed-layer  temperature  and  height,  9^ 
and  h,  respectively,  and  the  temperature  jump  When  there  is  no 

overshoot,  A0^  *  0,  and  the  mixed-layer  only  "encroaches"  on  the  elevated 
stable  layer. 


The  growth  of  the  convective  mixing  height  is  assumed  to  be  controlled 
by  the  bombardment  of  the  stable  lid  by  thermals  originating  at  the  surface. 
Computationally,  the  incremental  (hourly)  change  of  the  area  under  the 
temperature  profile  curve  (see  Figure  17)  is  proportional  to  the  hourly 
surface  heat  flux,  and  the  area  under  the  temperature  profile  9s  curve  is 
proportional  to  the  accumulated  surface  heat  flux.  Thus,  the  area  under  the 
temperature  profile  curve  at  time  t  is  given  by  the  formula: 


he  (h) 

s 


-f 


0  dz  a  (1  +  2A) 
s 


I 


t  <yx] 

pc. 


dr 


(47) 


2 

where  is  the  surface  heat  flux,  W/m  , 
p  Is  the  air  density, 

0^  is  the  specific  heat  of  air  at  constant  pressure,  and 
A  is  the  proportionality  constant  (ration  of  heat  flux  at  the  top  of 
boundary  layer  to  that  at  the  surface),  taken  as  0.2  after 
Deardor f f  (1980). 


Weil  and  Brower’s  formulation  of  the  mechanical  mixing  height,  after 
Tennekes  (1973)  and  Kato  and  Phillips  (1969),  is  similar  to  Eqn.  (47)  except 
that  the  right  side  of  the  equation  is  a  function  of  the  friction  velocity 
u,: 


ph  BT  t 

9  (h)  -  2  z0  dz  =  2  — 
s  J  s  g  J 

o  o 


u^(r)dx 


(48) 


where  Tq  is  the  surface  temperature  and  B  is  a  constant  (~  0.5). 
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For  each  hour,  the  higher  of  the  mechanical  or  convective  mixing  height 
is  chosen.  This  method  does  not  use  surface  temperature  explicitly  (unlike 
the  other  methods),  and  therefore  avoids  inconsistencies  between  the  on-site 
temperature  and  the  off-site  upper  sounding  temperature  at  122. 

The  interpolation  scheme  developed  by  Nieuwstadt  (1981)  for  the 
calculation  of  the  nocturnal  boundary  layer  height  (h)  is  used  in  SIGPRO: 

h  »  (L/3.8)  [-1+0+2. 28  u,/fL)1/2]  (49) 

where 


f  is  the  Coriolis  parameter,  =  2£2&in$; 

Q  is  the  angular  speed  of  rotation  of  the  earth,  =  7.292  x  10~5  rad  s' 
and 

0  is  the  latitude. 

1/2 

The  solution  for  h  given  by  this  formula  approaches  0.4(u»L/f)  for  small 
L,  as  suggested  by  Zilltinkevich  (1972),  and  approaches  0.  3u„/f  for  large  L 
(the  solution  for  neutral  conditions). 


2.4.4  Convective  Velocity  Scale 


In  the  convective  boundary  layer,  the  appropriate  velocity  scale  is  w#, 
which  can  be  computed  directly  from  its  definition  using  previously  defined 
estimates  of  and  h: 


= 


fg  Qh  h) 


T  p  c 


PJ 


1/3 


(50) 


If  water  vapor  buoyancy 
with  Qhv  and  T  replaced  with 


effects  are  important,  should  be  replaced 
Ty  in  Eqn.  (50),  as  described  in  Section  2.4.2. 
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3.  URBAN  DISPERSION  AND  STABILITY 


3. 1  Stability  Classification 

Most  dispersion  models  use  the  Pasquill  stability  classification  method 
(Pasquill,  1961)  or  a  variation  (e. g. ,  Turner,  1964)  to  define  atmospheric 
stability  and  dispersion  rates.  The  Pasquill  scheme  contains  six  discrete 
stability  classes  ranging  from  extremely  unstable  (class  A)  to  moderately 
stable  (class  F).  The  Pasquill  class  is  determined  based  on  wind  speed, 
solar  radiation,  and  cloudiness.  Solar  radiation  is  categorized  as  strong, 
moderate,  or  slight  (see  Table  5). 

The  Turner  (1964)  scheme  is  widely  used  in  current  EPA  dispersion 
models.  It  includes  a  modification  for  a  seventh  category  (extremely  stable 
-  class  7)  and  distinguishes  between  low  overcast  and  high  overcast 
conditions  before  assigning  the  neutral  stability  category.  The  relationship 
between  Turner's  classes  1-7  and  Pasquill’ s  classes  A-F  are  (Gifford,  1976; 
Hanna  et  al. ,  1982):  A«l,  B*2.  03,  D»4  and  5,  E=6,  and  F*7. 

Golder  (1972)  developed  a  nomogram  related  the  surface  roughness  length 
and  inverse  Monin-Obukhov  length  to  the  Pasquill  stability  classes.  The 
Golder  nomogram  is  shown  in  Figure  18.  Various  empirical  equations  have  been 
fit  to  Golder' s  graph  relating  L  to  zq  as  a  function  of  stability  class, 
including  equations  by  Shir  and  Shieh  (1974)  and  Irwin  (1979a) .  The  Irwin 
equation  was  based  on  a  power  law  relationship: 

L"1  =  a  (z  )b  (51) 

o 

where  the  coefficients  a  and  b  are  shown  in  Table  6. 

Lui  et  al.  (1976)  developed  an  empirical  expression  relating  Monin- 
Obukhov  length  to  surface  roughness  and  a  stability  parameter,  P  (Pasquill 
A  — ►  P  =  -3,  Pasquill  B  P  »  -2 . Pasquill  F  — »  P  =  +2),  such  that 
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Table  5 


The  Pasquill  (1961)  Stability  Categories 


10-m  wind 

speed  (m/s) 

Strong 

Insolation 

Moderate 

Slight 

-  Night 

Thinly  overcast  or 

>  4/8  low  cloud 

<  3/8  cloud 

<  2 

A 

A-B 

B 

F 

F 

2-3 

A-B 

B 

C 

E 

F 

3-5 

B 

B-C 

C 

D 

E 

5-6 

C 

C-D 

D 

D 

D 

>  6 

C 

D 

D 

D 

D 

Neutral  (D)  category  assumed  for  overcast  conditions,  day  or  night. 
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Figure  18.  Pasquill  stability  classes  as  a  function  of  inverse  Monin-Obukhov 
length  (1/L)  and  surface  roughness  length  (z^).  {Golder,  1972). 
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Table  6 


Power  Law  coefficients  Relating  L  to  zq  in 

the  Equation  (L  *  =  a  (z  )^) 

o 

(From  Irwin,  1979} 


Pasquill  Class  a 


b 


A 

B 

C 

D 

E 

F 


-0. 11360 
-0. 03849 
-0. 00807 
0.0 

0. 00807 
0. 03849 


-0. 1029 
-0. 1714 
-0.  3049 
0.0 

-0. 3049 
-0. 1714 
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L  « 


.  D  D3,  -(b.  -  b_ I P I  ♦  b_P2)' -1 

(a.  P  +  a_  P  )  z  1  2  3 

1  2  o 


(52  J 


where  the  empirical  constants  are  a^  =  0.004349,  a^  =  0.003724,  =  0.5034, 

b2  *  0.2310,  and  b3  =  0.0325. 

Bowling  (1985)  pointed  out  the  importance  of  snow  cover  in  determining 
atmospheric  stability  at  high  latitudes.  Also,  the  effect  of  latitude  on 
solar  eievation  angle  on  Pasquill  stability  was  determined  by  Andrews  and 
Hansen  (1988).  Hansen  (1990)  developed  a  fractional  stability  scheme 
combining  which  provides  greater  resolution  to  the  subjective  Pasquill 
classes  as  well  as  a  method  to  estimate  the  Monin-Obukhov  length.  This 
scheme  has  been  incorporated  into  the  FSCAT  computer  model  (Pena  and  Sutter, 
1990). 


Another  measure  of  stability  is  the  Kazanski -Monin  parameter,  p,  defined 
as  (Kazanski  and  Monin,  1960): 


8  *2  Qh 
f  p  cp  T  u2 


where  g 
k 

Qh 

f 

P 

c 

P 

T 

u* 


is  the  acceleration  due  to  gravity, 
is  the  von  Karman  constant, 
is  the  surface  sensible  heat  flux, 
is  the  Coriolis  parameter, 
is  the  air  density, 

is  the  specific  heat  of  air  at  constant  pressure, 
is  the  ambient  temperature,  and, 
is  the  surface  friction  velocity. 


(53) 


Smith  (1979)  developed  an  empirical  relationship  relating  p  to  Pasquill 
class,  P  for  unstable  and  neutral  cases  and  for  a  surface  roughness  length  of 
10  cm: 

3.  6 

P  =  - -  (54) 

1  -  0.  53(p/100 )  «•  4.9(p/100)  ♦  2.  0 (p/100 ) 
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Sutherland  et  al.  (1986)  defined  a  modified  version  of  the 
Kazanski-Monin  parameter,  p' ,  which  included  an  explicit  surface  roughness 
factor: 

g  k2  Qh 

p'  *  -a(z  )  - =-  (55) 

°  2fl  p  c  T  u, 

P 

where  a(z  )  =»  [a  ln(b/z  )]  A  and  a,  b  are  empirical  constants, 
o  o 

Sutherland  et  al.  (1986)  developed  a  simple  exponential  relationship 
relating  p'  to  the  Pasquill  class  P: 

P  =  A  eM'  (56) 

The  modified  Kazanski-Monin  parameter  eliminates  problems  with  p  at  low 
latitudes  associated  with  the  Coriolis  parameter  in  the  denominator  of  Eqn. 
(53),  and  is  more  general  than  Eqn.  (54)  in  that  it  applies  to  any 
stability  class  and  roughness  length. 

Sedefian  and  Bennett  (1980)  compared  the  results  of  several  schemes  for 
estimating  atmospheric  stability,  including  methods  based  on  the  standard 
deviation  of  horizontal  wind  fluctuations  (<r  method),  temperature  difference 

v7 

at  two  heights  (AT  method),  Richardson  number  and  bulk  Richardson  number 
methods,  and  Turner  (1964)  method.  Their  results  indicated  a  poor 
correlation  between  the  various  methods  on  an  hourly  basis,  and  they 
suggested  the  need  to  relate  plume  dispersion  coefficients  directly  to  the 
parameters  used  to  turbulence  parameters. 

Weil  and  Brower  (1984)  use  the  ratio  of  the  mean  wind  speed  to  the 
convective  velocity  scale  (i.e. ,  u/w.)  as  a  stability  parameter  rather  than 
the  Turner  class  during  convective  conditions.  Their  model,  the  Maryland 
Power  Plant  Siting  Program  (PPSP)  model  applies  to  tall  stacks  over  flat 
terrain.  They  use  rural  Briggs  dispersion  coefficients  for  elevated  sources 
in  rural  areas  with  the  ratio  u/w.  determining  the  effective  stability  class. 
During  nighttime  hours,  the  Turner  stability  class  scheme  is  used. •  Their 
stability  criteria  are  shown  in  Table  7.  In  urban  areas,  the  same  daytime 
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Table  7 

Stability  Scheme  for  Selecting  Dispersion  Parameters 
in  the  PPSP  Model  for  Rural  Areas 
(From  Weil  and  Brower,  1984) 


Day 


u/w.  criteria  Briggs  curve 


u/w.  <  3. 5  A 

3.  5  <  u/w.  <6.0  B 

6.  0  <  u/w.  <  13. 9  C 

13.9  <  u/w.  D 


Night  Turner  stability  class  Briggs  curve 


D 

E 

F 

G 


D 

E 

F 

F 
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stability  scheme  is  used  with  the  urban  Briggs  curves,  but  at  night,  urban 
Briggs  D  curves  are  applied  for  all  nighttime  Turner  stability  classes.  Some 
of  the  concepts  used  by  Weil  and  Brower  are  discussed  in  more  detail  in 
Section  3.2. 

Briggs  (1988)  also  concluded  that  the  ratio  u/w.  is  a  useful  measure  of 
dispersion  rates  in  moderately  to  very  unstable  conditions.  He  suggested 
that  convective  scaling  provides  the  best  framework  for  evaluating  the 
effects  of  urban  scale  surface  Inhomogeneities  on  dispersion  rates  in 
convective  conditions. 

Segal  et  al.  (1990)  examined  the  effects  of  surface  evaporation  rates  on 
Pasquill  classes  and  the  scaling  parameter  u/w..  Their  results  suggested 
that  increased  surface  moisture  has  a  significant  imprct  on  shifting 
stability  toward  more  neutral  conditions  as  measured  by  either  stability 
index.  High  surface  moisture  also  resulted  in  significantly  lower  daytime 
mixing  heights  (by  a  factor  of  1.5  to  2.0)  than  low  surface  moisture 
conditions.  These  results  are  attributed  to  a  shifting  of  available  energy 
away  from  sensible  heat  into  latent  heat  with  higher  values  of  surface 
moisture. 

The  EPA  Industrial  Source  Complex  (ISC)  model  accounts  for  differences 
in  surface  characteristics  on  dispersion  rates  by  employing  two  sets  of 
dispersion  curves:  the  Pasquill -Gifford  (PG)  curves  for  rural  areas  and  the 
McElroy-Pooler  (MP)  curves  for  urban  areas.  The  MP  curves  are  based  on  field 
experiments  conducted  in  the  St.  Louis  urban  area.  The  Briggs  curve  (Gifford, 
1976)  fit  to  the  St.  Louis  tracer  data  are  used  to  define  the  urban 
dispersion  coefficients  in  the  model.  The  Turner  (1964)  stability 
classification  scheme  is  used  in  ISC  with  either  PG  or  MP  dispersion 
coefficients.  In  urban  mode  with  the  MP  coefficients,  ISC  combines  stability 
classes  A-B  and  E-F,  resulting  in  four  separate  dispersion  classes.  In  rural 
mode,  six  stability  classes  (A-F)  are  used.  The  EPA  procedure  is  to  combine 
Turner  classes  6  and  7  into  class  F. 

One  of  the  difficulties  with  the  use  of  a  discrete  stability  class 
scheme  and  dispersion  curves  to  characterize  dispersion  is  the  discontinuous 
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nature  of  the  predicted  concentrations  as  a  function  of  stability  class  and 
surface  type.  In  addition,  because  the  PG  curves  were  derived  from 
near-ground. level  releases  over  relatively  flat  terrain  (z  -  3  cm)  and 
measured  only  over  short  downwind  distances  (i.e. ,  <  1  km),  their  application 
to  conditions  outside  this  range  of  conditions,  especially  to  plumes  emitted 
from  tall  stacks,  is  questionable  (Hanna  et  al. ,  1977). 

In  recent  years,  the  trend  in  dispersion  modeling  has  been  to 
parameterize  turbulence  (tnd  dispersion)  directly  in  terms  of  measured  or 
derived  meteorological  variables  using  boundary  layer  theory.  For  example, 
the  AMS  workshop  of  stability  classification  schemes  (Hanna  et  al. , 

1977)  suggested  the  use  of  the  standard  deviations  of  the  crosswind  component 
(c rv)  and  vertical  component  (cr^)  of  the  wind  to  estimate  the  plume  dispersion 
parameters  <r  and  o*^: 

<r  «  »  tf  (57 ) 

y  v  y 

o*  »  <r  t  f  (58) 

z  w  z 


where  t  is  time  for  the  pollutant  to  travel  to  the  receptor, 

<r  ,  <r  are  the  standard  deviations  of  the  horizontal  and  vertical 

y  z 

components  of  the  pollutant  concentration  distribution,  respectively, 

f  ,  f  are  nondlmensional  functions, 
y  z 


In  the  next  section,  methods  to  estimate  <r  ,  <t  , 

v  w 

boundary  layer  theory  are  described. 


f  ,  and  f  based  on 

y  z 


3.2  Similarity  Relations  for  Atmospheric  Turbulence  and  Dispersion 

Parameters 

A  number  of  different  scaling  laws  have  been  developed  to  apply  over 
various  regions  of  the  atmospheric  boundary  layer  or  during  particular 
stability  conditions.  The  basic  approach  is  to  represent  boundary  layer 
properties  in  terms  of  a  limited  number  of  characteristics  scaling 
parameters.  Holtslag  and  Nleuwstadt  (1986)  have  summarized  the  various 
regimes  and  the  scaling  parameters  which  apply  in  terms  of  two  diagrams 
(see  Figure  19  for  unstable  conditions  and  Figure  20  for  the  stable  boundary 
layer).  The  scaling  parameters  used  in  the  figures  are: 
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L  -  Monin-obukhov  length, 
h  -  mixing  height, 
z  -  height  above  the  ground, 
w8q  -  surface  height  flux, 
tq  -  surface  stress  (-  pu^),  and, 

A  -  local  length  scale  in  the  stable  boundary  layer 

Weil  (1985)  and  Briggs  (1985)  provide  reviews  on  the  use  of  similarity 
theory  in  diffusion  models.  In  the  convective  boundary  layer,  Weil  describes 
the  turbulence  characteristics  in  three  layers: 

(1)  Surface  layer:  z  a  0.  1  h  tr  ~  constant  with  height 

o*  increases  with  height 

(2)  Mixed  layer:  0.  i  h  <  z  <  0.8  h  or  ~  constant  with  height 

<r  ~  constant  with  height 

(3)  Entrainment  layer:  z  >  Q.8h  tr decreases  with  height 

<r^  decreases  with  height. 

(1977)  propose  the  following 

(59) 

(60) 

Hicks  (1985)  suggests  the  following  for  the  convective  mixed  layer  (0.1 
to  0.8  h). 


>  1/2 

3.  6 

2 

u. 

+  0.35 

2 

w* 

i 

(61  ) 

_  .  1/2 

1.2 

c. 

u. 

+  0.35 

C. 

w. 

(62) 
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In  the  neutral  boundary  layer,  Arya  (1984)  reports  monotonically 
decreasing  values  of  <r  and  <rw  throughout  the  mixed  layer.  Using  Blackadar 
and  Tennekes  (1968)  relationship  for  the  neutral  boundary  layer  height, 
Arya’ s  results  can  be  expressed  as: 


,  .  13e-0.9(z/h) 

v 


(63) 


(f 

w 


.  -  -0. 9(z/h) 

1.3  e 


(64) 


In  the  stable  boundary  layer,  Nieuwstadt  (1984)  finds  that  <ry  and  <r 
bear  constant  ratios  with  the  local  friction  velocity. 


Vu*l  *  Cv  (65) 

Vu-J  *  Cw  1661 

where  is  the  local  friction  velocity,  and, 

C  and  C  are  constants, 
v  w 

Hanna  et  al.  (1986)  suggest  that  C  *  1.6.  C  has  a  value  *  1.3 

v  w 

(Nieuwstadt,  1984).  The  local  friction  velocity,  u,j,,  can  be  expressed 
(Nieuwstadt,  1984)  as: 

3/4 

3  u.  (1  -  z/h)  (67) 

The  above  equations  represent  the  scaling  relations  appropriate  for  each 
layer  for  convective,  neutral,  and  stable  conditions.  However,  in  order  to 
provide  a  smooth  transition  from  one  vertical  layer  or  stability  regime  to 
another,  a  set  of  combined,  interpolation  formulas  have  been  developed  (Scire 
et  al,  ,  1990b).  These  equations  yield  the  proper  values  and  vertical 
variations  In  the  convective,  neutral,  and  stable  limits  while  providing  a 
mechanism  for  interpolating  the  results  for  intermediate  conditions  without 
physically  unrealistic  discontinuities.  The  following  equations,  based  on 
the  relations  presented  above,  apply  for  the  neutral-convective  boundary 
layer.  The  formulation  for  the  entrainment  layer  is  based  on  data  reported 
by  Caughey  (1981). 
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Surface  Layer:  z  a  0. 1  h 


(L  a  0) 


1/ 

£  4  +  0.35  w \ 

[  *n  '  2'9  "•  (-i/L)2/3] 


-0.  9(z/h) 


Mixed  Layer:  2  «  0.1  h  to  0.8  h  (L  a  0) 


rv  *  [  4  “•  \  *  0  35  “•  ] 
r„  -  [  1.15  i4  a2  *  0.35  w2  ] 


Entrainment  Layer:  z  >  0.8  h 


'*  -  [  4  U-  <  *  °':15  “•  ] 


(L  a  0) 


(68) 

(69) 

(70) 

(71) 

(72) 

(73) 


For  z  =  0.  8  h  to  1.2  h: 

r  2  2  2  ] 1/2 

<r  »  1.  15  u,  a4  +  a  ,  0.35  wt 

w  [  *  n  cl  •  J 

ac!  ~  [  0-5  +  (h-z)/(0. 4h) 


(74) 

(75) 

(76) 

(77) 


In  the  neutral-stable  boundary  layer,  the  'following  equations  are  used  t 
interpolate  vertical  profiles  of  <r^  and  as  a  function  of  stability.  They 
provide  the  proper  values  in  both  the  neutral  and  stable  limits. 


<r 

V 

=  u. 

[ {1-6  CS 

(z/L)  +1.8  an)/(l  +  z/L)J 

(L  >  0) 

(78) 

<r 

w 

=  1.3 

“•  [  (Cs 

(z/L)  +  an)/(l  +  z/L) j 

(L  >  0) 

(79) 

C 

C 

=  (1  - 

z/h)3/4 

(L  >  0) 

(80) 

In  order  to  provide  for  non-zero  plume  growth  rates  above  the  mixing 
height  and  to  prevent  numerical  problems  within  the  dispersion  calculations 
associated  with  near-zero  plume  dimensions,  minimum  <r  and  <r  values  can  be 

V  w 

assigned.  Hanna  et  ai.  ,  (1986)  suggest  that  an  appropriate  minimum  one-hour 
average  <ry  value  is  “  0.5  m/s. 


Equations  (68)  to  (77)  have  been  tested  with  the  original  data  providing 
the  basis  for  the  Panofsky  et  al.  (1977)  and  Hicks  (1985)  formulations.  The 
results,  summarized  in  Table  8,  indicate  that  the  modified  equations  compare 
well  with  the  original  equations  and  the  observational  data.  The  modified 
equations  have  the  advantage  of  a  smooth  and  continuous  transition  to  the 
neutral  results  of  Arya  (1984). 

Irwin  (1983)  has  evaluated  several  schemes  for  determining  the  f  and  f 

y  2 

functions,  including  equations  suggested  by  Draxler  (1976),  Pasquill  (1976), 
and  Cramer  (1976).  His  conclusion  was  that  the  parameterization  suggested  by 
Draxler  (1976)  performed  best  overall. 

fy  =  £  1  +  0.9  (t/1000)1/2j  (81) 

1+0.9  (t/500) 1/2 
1  +  0.945  (t/100)0-806  J 

These  relations  provide  a  mechanism  for  determining  atmospheric 
turbulence  and  dispersion  rates  in  a  continuous  manner  as  a  function  of 
stability  and  plume  height. 


(L  < 

0) 

(82) 

(L  > 

0) 

(83) 
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Table  8 

Comparison  of  Panofsky  et  ai.  (1977)  and  Hicks  (1985) 
Formulations  of  o*v  and  <ry  with  Equations  (68)-(77) 
(From  Scire  et  al.  ,  1990) 


Panofsky  et  al.  data 

Observed  <r  vs. 

V 

Observed  <r  vs. 

V 

Panofsky  <r  vs. 

Panofsky 

Eqns.  ( 68 ) — ( 77 ) 

Eqns.  (68) -(77) 

Averages 

(1.14.  1.20) 

(1.14,  1.21) 

(1.20,  1.21) 

Corr.  Coef. 

.  81 

.  84 

.  992 

Average  Bias 

.07 

.  07 

.00 

Average  Abs.  Error 

.  10 

.09 

.02 

RMSE 

.  13 

.  12 

.02 

Hicks  1985  data 

Observed  <r  vs. 

V 

Observed  <r  vs. 

V 

Hicks  <r  vs. 

V 

Hicks 

Eqns.  (68) -(77) 

Eqns.  (68) -(77) 

Averages 

(1.17,  1.12) 

(1.17,  1.06) 

(1.12,  1.06) 

Corr.  Coef. 

.79 

.77 

.  998 

Average  Bias 

-.05 

-.  11 

.06 

Average  Abs.  Error 

.20 

.23 

.  06 

RMSE 

.27 

.  30 

.08 

Hicks  1985  data 

Observed  <r  vs. 
w 

Hicks 

Observed  <r  vs. 
w 

Eqns.  (68)- (77) 

Hicks  <r  vs. 
w 

Eqns.  ( 68 )  —  ( 77 ) 

Averages 

(0.98,  1.01) 

(0.98,  0.98) 

(1.01,  0.98) 

Corr.  Coef. 

.  91 

.  91 

.  998 

Average  Bias 

.03 

.  00 

03 

Average  Abs.  Error 

.  12 

.  11 

.  03 

RMSE 

.  15 

.  14 

.  04 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  Phase  I  effort,  a  detailed  methodology  appropriate  for 
estimating  stability  and  dispersion  parameters  in  an  urban  environment  has 
been  described.  The  main  advantages  of  the  approach  include  the  ability  to 
provide  continuous  estimates  of  stability  and  dispersion  as  a  function  of 
stability,  plume  height,  and  surface  characteristics,  the  requirement  for 
only  routinely-available  meteorological  inputs  and  land  use  data,  and 
computational  efficiency  (i.e. ,  it  is  designed  to  run  on  a  PC). 

The  major  components  of  the  model  are  a  new  soil  moisture  model  capable 
of  predicting  surface  wetness  factors  over  both  rural  and  urban  land  use 
types,  a  meteorological  preprocessor  which  computes  heat  fluxes  and  boundary 
layer  scaling  parameters,  and  a  set  of  similarity  relationships  to  estimate 
atmospheric  stability  and  plume  dispersion  rates.  The  various  components  of 
the  model  have  been  coded  and  subjected  to  module  testing.  Among  the 
conclusions  and  recommendations  of  the  Phase  I  literature  review  and 
development  effort  are: 

(1)  The  surface  moisture  parameter  is  a  critical  input  to  the  energy 
balance  scheme  of  Holtslag  and  vap  Ulden  (1983).  It  is  also  one 
of  the  most  difficult  to  specify  because  it  exhibits  considerable 
spatial  and  time  variability,  reflecting  differences  In  land  use 
and  the  variability  of  precipitation. 

(2)  A  simple  two-capacitor  soil  moisture  model  has  shown  the  ability  to 
reproduce  the  behavior  of  more  complex  numerical  models  during 
simulated  drying  out  periods.  Sample  simulations  of  a  three  month 
period  for  both  rural  and  urban  areas  produce  qualitatively 
reasonable  results.  Further  testing  and  enhancement  of  the 
capacitance  model  with  the  Indianapolis  urban  data  base  is 
recommended.  The  testing  should  focus  on  two  issue:  1)  Do  observed 
heat  fluxes  show  the  expected  behavior  during  periods  with  no 
precipitation?  (i.e.,  a  shift  from  latent  to  sensible  heat  flux 
with  time);  2)  Is  the  capacitance  model  able  to  more  accurately 
predict  energy  fluxes  than  the  use  of  constant,  time-averaged 
values  of  the  soil  moisture  parameter? 
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(3) 


The  soil  moisture  model  contains  several  parameters  which  need  to 
be  determined  empirically.  These  include  the  capacitor  time 
constants  (t)  and  minimum/maximum  capacitor  change  (q).  Additional 
work  needs  to  be  done  to  (1)  evaluate  if  the  number  of  parameters 
can  be  reduced  and  (2)  determine  recommended  values  for  the 
parameters  based  on  routinely  available  information  (e.g.  ,  land  use 
types ) . 

(4)  A  method  for  including  the  buoyancy  effects  of  water  vapor  on 
boundary  layer  parameters  has  been  described.  Further  testing  is 
recommended  to  quantify  the  magnitude  of  these  effects  and  the 
conditions  under  which  they  are  likely  to  be  important. 

(5)  A  well-known  surface  energy  balance  method,  adapted  for  urban 
conditions,  but  also  applicable  to  suburban  and  rural  land  surfaces 
as  well,  has  been  demonstrated.  The  method  is  well  suited  for 
further  study  and  development.  It  requires  only  routinely 
available  data  inputs  and  produces  the  boundary  layer  parameters 
required  by  state-of-the-science  scaling  techniques  to  estimate 
stability  and  plume  dispersion  rates.  It  is  recommended  that  the 
performance  of  the  coupled  capacitance-meteorological  model  be 
assessed  in  future  Phase  II  work.  The  Indianapolis  and  St.  Louis 
data  sets  are  suitable  candidates  for  this  evaluation  effort. 

(6)  A  number  of  schemes  to  estimate  Pasquill  stability  classes  have  been 
reviewed.  The  use  of  scaling  techniques  is  recommended  to  provide 

a  continuous  estimate  of  dispersion  and  stability  parameters 
consistent  with  current  boundary  layer  similarity  theory. 

(7)  A  series  of  scaling  equations  describing  the  turbulence  parameters, 

< and  in  various  scaling  regimes  (height  and  stability 
regimes)  have  been  presented.  Two  sets  of  interpolation  formulas 
which  provide  for  a  smooth  and  continuous  variation  of  the 
turbulence  parameters  from  one  regime  to  another  have  been 
described.  The  first  set  deals  with  the  unstable  to  neutral 
regime,  and  the  second  handles  the  stable  to  neutral  transition. 
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Further  refinement  and  testing  of  these  formulas  is  recommended  in 
a  Phase  II  effort. 

(8)  Once  the  three  major  modules  of  the  model  have  been  tested  further, 
it  is  recommended  that  the  components  be  integrated  into  a  single 
system  and  tested  as  a  unit  in  Phase  II  work.  An  user-friendly 
front-end  input  module  should  be  included  to  facilitate  its  use. 

The  predicted  dispersion  parameters  are  suitable  for  coupling 
directly  to  a  dispersion  model. 

This  Phase  I  study  has  demonstrated  the  feasibility  of  a  pc-based 
modeling  system  for  predicting  urban  (and  rural)  stability  and  atmospheric 
turbulence  parameters  in  a  continuous  manner  based  on  the  use  of  only 
routinely-available  input  data.  The  actual  development  of  an  integrated 
system  and  further  testing  of  the  system  components  and  the  completed  unit 
are  recommended  as  a  potential  Phase  II  development  project. 
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